Siesta

Merge lp:~nickpapior/siesta/low-D-stress into lp:siesta/4.1

  1. low-D-stress
  2. Merge into rel-4.1
Proposed by Nick Papior
Status: Needs review
Proposed branch: lp:~nickpapior/siesta/low-D-stress
Merge into: lp:siesta/4.1
Diff against target: 2798 lines (+2016/-685) (has conflicts)
12 files modified
Src/Makefile (+5/-3)
Src/write_subs.F90 (+884/-681)
Tests/Makefile (+1/-1)
Tests/Reference/graphene_non_perp.out (+521/-0)
Tests/Reference/graphene_perp.out (+521/-0)
Tests/graphene_non_perp/graphene_non_perp.fdf (+33/-0)
Tests/graphene_non_perp/graphene_non_perp.pseudos (+1/-0)
Tests/graphene_non_perp/makefile (+6/-0)
Tests/graphene_perp/graphene_perp.fdf (+33/-0)
Tests/graphene_perp/graphene_perp.pseudos (+1/-0)
Tests/graphene_perp/makefile (+6/-0)
version.info (+4/-0)
Text conflict in version.info
To merge this branch: bzr merge lp:~nickpapior/siesta/low-D-stress
Reviewer Review Type Date Requested Status
Alberto Garcia Pending
Review via email: mp+363189@code.launchpad.net

Commit message

Added 1D and 2D stress tensor calculations

The problem with low-D materials is that stress is "per volume",
but since vacuum is independent in a slab calculation the stress will
be variable depending on the amount of vacuum.

What one *really* should do is to use the volume of the slab along the
vacuum direction so it may be fixed depending on the extend of the slab.
This poses additional problems since how do you decide the slab thickness?
- vdW bond-lengths
- orbital radii
- or 3rd?

This amends output to do the simplest thing but forces users to do post-processing.
The vacuum region is divided out and thus change the units depending on the dimensionality.
- 2D == eV/Ang^2
- 1D == eV/Ang

Then the user can them-selves determine the actual height/area of the slab/chain.

Note that the regular stress tensors are still printed out, but now an additional
block is added.

Description of the change

In addition to the above message I would like to discuss the following which came up during the coding:

1) I am currently using nsc to determine the periodic directions, not ideal but I put in this remark in the code:
    ! TODO, probably we need to use shaper...
    ! However, the shaper returns vectors which are
    ! not parallel to the actual lattice vectors and hence
    ! it becomes difficult to utilize them
    ! What to do for Gamma calculations?

I tried using shaper, to no avail as noted above...

2) Regarding forces in the write_forces. Currently ghost atoms are also considered in sum(fa) and res(fa), don't you agree we need to change that to not include ghost atoms (at least in the residual?) Or?

To post a comment you must log in.
Revision history for this message
Alberto Garcia (albertog) wrote :

A couple of preliminary questions:

1. Are the resulting 1D and 2D stresses vacuum independent now?
2. At least for surfaces, a physically meaningful concept is that of the
'surface stress', for which a computational recipe exists. I do not think
that your 2D stress will map to that in the case of a slab meant to
represent a surface. The concept of surface stress does not apply to
one-atomic-layer materials, however.

On Thu, Feb 14, 2019, 11:05 Nick Papior <<email address hidden> wrote:

> Nick Papior has proposed merging lp:~nickpapior/siesta/low-D-stress into
> lp:siesta/4.1.
>
> Commit message:
> Added 1D and 2D stress tensor calculations
>
> The problem with low-D materials is that stress is "per volume",
> but since vacuum is independent in a slab calculation the stress will
> be variable depending on the amount of vacuum.
>
> What one *really* should do is to use the volume of the slab along the
> vacuum direction so it may be fixed depending on the extend of the slab.
> This poses additional problems since how do you decide the slab thickness?
> - vdW bond-lengths
> - orbital radii
> - or 3rd?
>
> This amends output to do the simplest thing but forces users to do
> post-processing.
> The vacuum region is divided out and thus change the units depending on
> the dimensionality.
> - 2D == eV/Ang^2
> - 1D == eV/Ang
>
> Then the user can them-selves determine the actual height/area of the
> slab/chain.
>
> Note that the regular stress tensors are still printed out, but now an
> additional
> block is added.
>
> Requested reviews:
> Alberto Garcia (albertog)
>
> For more details, see:
> https://code.launchpad.net/~nickpapior/siesta/low-D-stress/+merge/363189
>
> In addition to the above message I would like to discuss the following
> which came up during the coding:
>
> 1) I am currently using nsc to determine the periodic directions, not
> ideal but I put in this remark in the code:
> ! TODO, probably we need to use shaper...
> ! However, the shaper returns vectors which are
> ! not parallel to the actual lattice vectors and hence
> ! it becomes difficult to utilize them
> ! What to do for Gamma calculations?
>
> I tried using shaper, to no avail as noted above...
>
>
> 2) Regarding forces in the write_forces. Currently ghost atoms are also
> considered in sum(fa) and res(fa), don't you agree we need to change that
> to not include ghost atoms (at least in the residual?) Or?
> --
> You are requested to review the proposed merge of
> lp:~nickpapior/siesta/low-D-stress into lp:siesta/4.1.
>

Reply
Revision history for this message
Nick Papior (nickpapior) wrote :
Download full text (4.4 KiB)

1)

For 2D graphene I find the following dependence on the 2D stress tensor (directory is vacuum in Ang)

100/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
100/RUN.out-siesta: -0.104105 -0.000000 0.000000
100/RUN.out-siesta: -0.000000 -0.104441 0.000000
--
110/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
110/RUN.out-siesta: -0.104470 -0.000000 0.000000
110/RUN.out-siesta: -0.000000 -0.104806 0.000000
--
120/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
120/RUN.out-siesta: -0.104281 -0.000000 0.000000
120/RUN.out-siesta: 0.000000 -0.104617 0.000000
--
130/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
130/RUN.out-siesta: -0.104602 -0.000000 0.000000
130/RUN.out-siesta: -0.000000 -0.104938 0.000000
--
20/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
20/RUN.out-siesta: -0.104098 -0.000000 0.000000
20/RUN.out-siesta: 0.000000 -0.104434 0.000000
--
30/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
30/RUN.out-siesta: -0.105210 -0.000000 0.000000
30/RUN.out-siesta: -0.000000 -0.105546 0.000000
--
40/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
40/RUN.out-siesta: -0.104277 0.000000 0.000000
40/RUN.out-siesta: 0.000000 -0.104613 0.000000
--
50/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
50/RUN.out-siesta: -0.105209 0.000000 0.000000
50/RUN.out-siesta: 0.000000 -0.105545 0.000000
--
60/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
60/RUN.out-siesta: -0.104109 -0.000000 0.000000
60/RUN.out-siesta: -0.000000 -0.104445 0.000000
--
70/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
70/RUN.out-siesta: -0.104656 0.000000 0.000000
70/RUN.out-siesta: 0.000000 -0.104992 0.000000
--
80/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
80/RUN.out-siesta: -0.104286 -0.000000 0.000000
80/RUN.out-siesta: 0.000000 -0.104622 0.000000
--
90/RUN.out:siesta: 2D stress tensor (static) (eV/Ang**2):
90/RUN.out-siesta: -0.105208 0.000000 0.000000
90/RUN.out-siesta: 0.000000 -0.105544 0.000000

For a 1D carbon chain with directory being the vacuum along the two perpendicular directions
10/RUN.out:siesta: 1D stress tensor (static) (eV/Ang):
10/RUN.out-siesta: 0.000000 0.000000 0.000000
10/RUN.out-siesta: 0.000000 0.000000 0.000000
10/RUN.out-siesta: 0.000000 0.000000 7.196726
--
20/RUN.out:siesta: 1D stress tensor (static) (eV/Ang):
20/RUN.out-siesta: 0.000000 0.000000 0.000000
20/RUN.out-siesta: 0.000000 0.000000 0.000000
20/RUN.out-siesta: 0.000000 0.000000 7.198956
--
30/RUN.out:siesta: 1D stress tensor (static) (eV/Ang):
30/RUN.out-siesta: 0.000000 0.000000 0.000000
30/RUN.out-siesta: 0.000000 0.000000 0.000000
30/RUN...

Read more...

Reply
lp:~nickpapior/siesta/low-D-stress updated
1071. By Nick Papior

Added two examples where this should be printed out

Currently I have not added the output since it may change.

1072. By Nick Papior

Finalized calculation of low-D stress tensors using the shaper routine

We now calculate the vacuum regions based on the shaper routine.
This results in very similar output for skewed lattice vectors.
See e.g. the two new tests.

  graphene_perp
  graphene_non_perp

Revision history for this message
Alberto Garcia (albertog) wrote :

For variable-cell relaxations, is it possible to use as 'target-stress' a
similarly renormalized value?

Reply
Revision history for this message
Nick Papior (nickpapior) wrote :

Currently not.

What we could do is to read in a TargetStress1D/2D and use this (renormalized) to TargetStress?

In that way we don't have to change any of the other routines? However, I would think that this is more a special case since one can always *pre-calculate* this from outside in the same way we do inside?

Reply

Unmerged revisions

1072. By Nick Papior

Finalized calculation of low-D stress tensors using the shaper routine

We now calculate the vacuum regions based on the shaper routine.
This results in very similar output for skewed lattice vectors.
See e.g. the two new tests.

  graphene_perp
  graphene_non_perp

1071. By Nick Papior

Added two examples where this should be printed out

Currently I have not added the output since it may change.

1070. By Nick Papior

Added 1D and 2D stress tensor calculations

The problem with low-D materials is that stress is "per volume",
but since vacuum is independent in a slab calculation the stress will
be variable depending on the amount of vacuum.

What one *really* should do is to use the volume of the slab along the
vacuum direction so it may be fixed depending on the extend of the slab.
This poses additional problems since how do you decide the slab thickness?
- vdW bond-lengths
- orbital radii
- or 3rd?

This amends output to do the simplest thing but forces users to do post-processing.
The vacuum region is divided out and thus change the units depending on the dimensionality.
- 2D == eV/Ang^2
- 1D == eV/Ang

Then the user can them-selves determine the actual height/area of the slab/chain.

Note that the regular stress tensors are still printed out, but now an additional
block is added.

TODO: these stress's should probably be the ones used for relaxations and var-cell
      calculations?

1069. By Nick Papior

Cleaned up write_subs.F90 and changed implicit loops to explicit notations

This merely makes the code more readable and makes it ready for the low-D changes.

Fixed all use statements in write_subs.
Moved siesta_write_forces to a separate module, along the lines of the pressure,
energetics, etc.

1068. By Nick Papior

Moved write_subs.F to write_subs.F90 for easier maintenance

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'Src/Makefile'
2--- Src/Makefile 2019-02-13 17:39:22 +0000
3+++ Src/Makefile 2019-02-28 08:28:01 +0000
4@@ -1341,9 +1341,10 @@
5 write_md_record.o: md_out.o parallel.o siesta_geom.o siesta_options.o units.o
6 write_orb_indx.o: atmfuncs.o cellsubs.o files.o precision.o
7 write_raw_efs.o: atmfuncs.o atomlist.o precision.o siesta_geom.o
8-write_subs.o: alloc.o atomlist.o m_energies.o m_forces.o m_iostruct.o m_spin.o
9-write_subs.o: m_steps.o m_stress.o m_ts_global_vars.o parallel.o precision.o
10-write_subs.o: siesta_cml.o siesta_geom.o siesta_options.o units.o zmatrix.o
11+write_subs.o: atomlist.o intrinsic_missing.o m_energies.o m_forces.o
12+write_subs.o: m_iostruct.o m_spin.o m_steps.o m_stress.o m_ts_global_vars.o
13+write_subs.o: parallel.o precision.o siesta_cml.o siesta_geom.o
14+write_subs.o: siesta_options.o units.o zmatrix.o
15 writewave.o: alloc.o atmfuncs.o atomlist.o densematrix.o diag.o diag_option.o
16 writewave.o: get_kpoints_scale.o kpoint_grid.o m_spin.o parallel.o
17 writewave.o: parallelsubs.o precision.o siesta_geom.o sys.o units.o
18@@ -1507,5 +1508,6 @@
19 trialorbitalclass.o: m_trialorbitalclass.o
20 version_info.o: version.o
21 write_subs_energies.o: write_subs.o
22+write_subs_forces.o: write_subs.o
23 write_subs_positions.o: write_subs.o
24 write_subs_pressure.o: write_subs.o
25
26=== renamed file 'Src/write_subs.F' => 'Src/write_subs.F90'
27--- Src/write_subs.F 2017-11-23 14:17:27 +0000
28+++ Src/write_subs.F90 2019-02-28 08:28:01 +0000
29@@ -11,698 +11,901 @@
30 ! The submodules appear first in this file, and module
31 ! write_subs is at the end.
32 !
33- Module write_subs_positions
34-
35- use m_steps, only: final
36-
37- private
38- public :: siesta_write_positions
39- CONTAINS
40-
41- subroutine siesta_write_positions(moved)
42- use siesta_geom
43- use atomlist, only: elem, iza
44- use siesta_cml
45- use units, only: Ang
46- use zmatrix, only: lUseZmatrix, write_canonical_ucell_and_Zmatrix
47- use m_iostruct, only: write_struct
48- use siesta_options, only: idyn
49-
50- implicit none
51-
52-! Whether the structure has already changed in response to forces and stresses
53- logical, intent(in) :: moved
54-
55-!
56-! Write out structural information in "crystallography" format,
57-! and in "canonical" (see Zmatrix module) zmatrix format.
58-
59-
60- call write_struct( ucell, na_u, isa, iza, xa, moved)
61- if (lUseZmatrix .and. (idyn .eq. 0)) then
62- if (moved) then
63- call write_canonical_ucell_and_Zmatrix
64- $ (filename="NEXT_ITER.UCELL.ZMATRIX")
65- else
66- call write_canonical_ucell_and_Zmatrix
67- $ (filename="OUT.UCELL.ZMATRIX")
68- endif
69- endif
70-
71- ! Note that this information should be written only at the
72- ! time of processing the current geometry (in state_init)
73-
74- if ((.not. moved) .and. cml_p) then
75- call cmlAddMolecule(xf=mainXML, natoms=na_u, elements=elem,
76- . atomRefs=cisa, coords=xa/Ang)
77- call cmlAddLattice(xf=mainXML, cell=ucell/Ang,
78- . units='siestaUnits:Ang', dictref='siesta:ucell')
79- endif
80-
81- end subroutine siesta_write_positions
82-
83- END MODULE write_subs_positions
84- Module write_subs_pressure
85-
86- use m_steps, only: final
87-
88- private
89- public :: siesta_write_stress_pressure
90-
91- CONTAINS
92-
93- subroutine siesta_write_stress_pressure()
94-
95- use parallel, only: IOnode
96- use precision
97- USE siesta_options
98- use siesta_geom
99- use atomlist, only: iza
100- use m_iostruct, only: write_struct
101- use siesta_cml
102- use units
103- use m_spin
104- use m_energies, only: FreeE
105- use m_stress, only: stress, kin_stress, mstress, cstress, tstress
106-
107- implicit none
108-
109- integer :: jx, ix
110-
111- real(dp):: Pmol ! Molecular pressure (discounting Virial term)
112- real(dp):: Psol ! Pressure of "solid"
113- real(dp):: Press ! Pressure
114- real(dp):: ps(3,3) ! Auxiliary array
115-
116-! Stress tensor and pressure:
117-
118- if (.not.final) then
119-!
120-! Write Voigt components of total stress tensor
121-!
122- ps = stress + kin_stress
123- write(6,'(/,a,6f12.2)')
124- . 'Stress-tensor-Voigt (kbar):',
125- . (ps(jx,jx)/kbar,jx=1,3),
126- $ ps(1,2)/kbar,
127- $ ps(2,3)/kbar,
128- $ ps(1,3)/kbar
129- Press = - ((ps(1,1) + ps(2,2) + ps(3,3))/3.0_dp)
130- write(6,"(a,f14.4)") "(Free)E + p*V (eV/cell)",
131- $ (FreeE + Press*volume_of_some_cell)/eV
132-
133- if (RemoveIntraMolecularPressure) then
134- ps = mstress + kin_stress
135- write(6,'(/,a,6f12.2)')
136- . 'Inter-Molecular-Stress-Voigt (kbar):',
137- . (ps(jx,jx)/kbar,jx=1,3),
138- $ ps(1,2)/kbar,
139- $ ps(2,3)/kbar,
140- $ ps(1,3)/kbar
141- Press = - ((ps(1,1) + ps(2,2) + ps(3,3))/3.0_dp)
142- write(6,"(a,f14.4)")
143- $ "(Free)E + p_inter_molec * V (eV/cell)",
144- $ (FreeE + Press*volume_of_some_cell)/eV
145- endif
146-!
147-! This use of the volume is OK, as it is called from state_analysis,
148-! before possibly changing the cell.
149-! Write "target enthalpy" (E + pV, where p is the *target* pressure)
150- write(6,"(a,f14.4)") "Target enthalpy (eV/cell)",
151- $ (FreeE + tp*volume_of_some_cell)/eV
152-
153- ! Write stress to CML file always
154- if (cml_p) then
155- call cmlAddProperty(xf=mainXML, value=stress*Ang**3,
156- . dictref='siesta:stress', title='Stress',
157- . units='siestaUnits:evpa3')
158- endif !cml_p
159-
160- ! Output depends on dynamics option
161- select case (idyn)
162- case(0:5,8,9)
163-
164- if (idyn==0 .and. (.not.varcel)) then
165- continue
166- else
167- write(6,'(/,a,3(/,a,3f12.6))')
168- . 'siesta: Stress tensor (static) (eV/Ang**3):',
169- . (' ',(stress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
170- Psol = - ((stress(1,1) + stress(2,2) + stress(3,3))/3.0_dp)
171- write(6,'(/,a,f20.8,a)')
172- . 'siesta: Pressure (static):', Psol/kBar, ' kBar'
173-! write(6,'(/,a,3(/,a,3f12.6))')
174-! . 'siesta: Stress tensor (static-constrained) (eV/Ang**3):',
175-! . (' ',(cstress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
176-! write(6,'(a,f20.8,a)')
177-! . 'siesta: Pressure (static-constrained):',
178-! . - ((cstress(1,1)+cstress(2,2)+cstress(3,3))/3.0_dp)/kBar,
179-! . ' kBar'
180- if (cml_p) then
181- ! stress written above
182- call cmlAddProperty(xf=mainXML, value=Psol,
183- . dictref='siesta:psol', title='Pressure (Static)',
184- . units='siestaUnits:kBar')
185- endif !cml_p
186- write(6,'(/,a,3(/,a,3f12.6))')
187- . 'siesta: Stress tensor (total) (eV/Ang**3):',
188- . (' ',(tstress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
189- Psol = - ((tstress(1,1)+tstress(2,2) +tstress(3,3))/3.0_dp)
190- write(6,'(/,a,f20.8,a)')
191- . 'siesta: Pressure (total):', Psol/kBar, ' kBar'
192- if (cml_p) then
193- call cmlAddProperty(xf=mainXML, value=tstress*Ang**3,
194- . dictref='siesta:tstress', title='Total Stress',
195- . units='siestaUnits:evpa3')
196- call cmlAddProperty(xf=mainXML, value=Psol,
197- . dictref='siesta:tpsol', title='Pressure (Total)',
198- . units='siestaUnits:kBar')
199- endif !cml_p
200-
201- if (RemoveIntraMolecularPressure) then
202- ps = mstress
203- write(6,'(/,a,3(/,a,3f12.6))')
204- . 'siesta: Stress tensor (nonmol) (eV/Ang**3):',
205- . (' ',(ps(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
206- Psol = - ((ps(1,1) + ps(2,2) + ps(3,3))/3.0_dp)
207- write(6,'(/,a,f20.8,a)')
208- . 'siesta: Pressure (nonmol):', Psol/kBar, ' kBar'
209- if (cml_p) then
210- call cmlAddProperty(xf=mainXML, value=ps*Ang**3,
211- . dictref='siesta:mstress',
212- . title='Stress tensor (normal)',
213- . units='siestaUnits:evpa3')
214- call cmlAddProperty(xf=mainXML, value=Psol,
215- . dictref='siesta:pmol', title='Pressure (Nonmol)',
216- . units='siestaUnits:kBar')
217- endif !cml_p
218-!
219- ps = mstress + kin_stress
220- write(6,'(/,a,3(/,a,3f12.6))')
221- . 'siesta: Stress tensor (nonmol+kin) (eV/Ang**3):',
222- . (' ',(ps(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
223- Psol = - ((ps(1,1)+ps(2,2) +ps(3,3))/3.0_dp)
224- write(6,'(/,a,f20.8,a)')
225- . 'siesta: Pressure (nonmol+kin):', Psol/kBar, ' kBar'
226- if (cml_p) then
227- call cmlAddProperty(xf=mainXML, value=ps*Ang**3,
228- . dictref='siesta:tmstress',
229- . title='Stress tensor (nonmol+kin)',
230- . units='siestaUnits:evpa3')
231- call cmlAddProperty(xf=mainXML, value=Psol,
232- . dictref='siesta:tpmol', title='Pressure (Nonmol+Kin)',
233- . units='siestaUnits:kBar')
234- endif !cml_p
235- endif ! Remove intramolecular pressure
236-
237- endif !varcel
238-
239- end select !idyn
240-
241- else !final
242-
243- ! AG: Possible BUG
244- ! The volume here refers to the "old" cell, and
245- ! might be out of date if the cell has changed.
246-
247-! Print stress tensor unconditionally
248- write(6,'(/,a,3(/,a,3f12.6))')
249- . 'siesta: Stress tensor (static) (eV/Ang**3):',
250- . ('siesta: ',(stress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
251- if (cml_p) then
252- call cmlAddProperty(xf=mainXML, value=stress*Ang**3/eV,
253- . dictref='siesta:stress', units='siestaUnits:eV_Ang__3')
254- endif !cml_p
255-
256-! Print constrained stress tensor if different from unconstrained
257- if (Any(cstress /= stress )) then
258- write(6,'(/,a,3(/,a,3f12.6))')
259- . 'siesta: Constrained stress tensor (static) (eV/Ang**3):',
260- . ('siesta: ',(cstress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
261- if (cml_p) then
262- call cmlAddProperty(xf=mainXML, value=cstress*Ang**3/eV,
263- . dictref='siesta:cstress',
264- . units='siestaUnits:eV_Ang__3')
265- endif !cml_p
266- endif
267-
268-! Find pressure
269-
270- Psol = - (( stress(1,1) + stress(2,2) + stress(3,3) )/3.0_dp)
271- Pmol = - (( mstress(1,1) + mstress(2,2) + mstress(3,3) )/3.0_dp)
272-
273- write(6,'(/,a,f18.6,a)')
274- . 'siesta: Cell volume =',
275- $ volume_of_some_cell/Ang**3, ' Ang**3'
276- write(6,'(/,a,/,a,2a20,a,3(/,a,2f20.8,a))')
277- . 'siesta: Pressure (static):',
278- . 'siesta: ','Solid', 'Molecule', ' Units',
279- . 'siesta: ', Psol, Pmol, ' Ry/Bohr**3',
280- . 'siesta: ', Psol*Ang**3/eV, Pmol*Ang**3/eV, ' eV/Ang**3',
281- . 'siesta: ', Psol/kBar, Pmol/kBar, ' kBar'
282- if (cml_p) then
283- call cmlStartPropertyList(mainXML, title='Final Pressure')
284- call cmlAddProperty(xf=mainXML,
285- $ value=volume_of_some_cell/Ang**3,
286- . title='cell volume', dictref='siesta:cellvol',
287- . units='siestaUnits:Ang__3')
288- call cmlAddProperty(xf=mainXML, value=Psol/kBar,
289- . title='Pressure of Solid', dictref='siesta:pressSol',
290- . units='siestaUnits:kbar')
291- call cmlAddProperty(xf=mainXML, value=Pmol/kBar,
292- . title='Pressure of Molecule', dictref='siesta:pressMol',
293- . units='siestaUnits:kbar')
294- call cmlEndPropertyList(mainXML)
295- endif !cml_p
296-
297- endif !final for stress & pressure
298-
299- end subroutine siesta_write_stress_pressure
300-
301-
302- End Module write_subs_pressure
303-
304- Module write_subs_energies
305-
306- use m_steps, only: final, istp
307-
308- private
309- public :: siesta_write_energies
310- CONTAINS
311-
312- subroutine siesta_write_energies( iscf, dDmax, dHmax )
313- USE siesta_options
314- use siesta_cml
315- use units
316- use m_energies
317- use m_spin
318- use m_ts_global_vars, only: TSinit, TSrun
319- implicit none
320-
321- integer :: iscf
322- real(dp), intent(in) :: dDmax ! Max. change in dens. matrix elem.
323- real(dp), intent(in) :: dHmax ! Max. change in H elements
324-
325- character(len=64) :: fmt
326- character(len=6) :: scf_name
327-
328- logical :: first_scf_step
329- integer :: i
330-
331- scf_name = 'scf'
332- if ( TSrun ) then
333- scf_name = 'ts-scf'
334- end if
335-
336- first_scf_step = (iscf == 1)
337- ! Only print out full decomposition at very beginning and end.
338- if ((istp==1.and.first_scf_step).or.final) then
339- write(6,'(/,a,/,(a,f17.6))')
340- . 'siesta: Program''s energy decomposition (eV):',
341- . 'siesta: Ebs =', Ebs/eV,
342- . 'siesta: Eions =', Eions/eV,
343- . 'siesta: Ena =', Ena/eV,
344- . 'siesta: Ekin =', Ekin/eV,
345- . 'siesta: Enl =', Enl/eV,
346- . 'siesta: Eso =', Eso/eV,
347- . 'siesta: Eldau =', Eldau/eV,
348- . 'siesta: DEna =', DEna/eV,
349- . 'siesta: DUscf =', DUscf/eV,
350- . 'siesta: DUext =', DUext/eV,
351- . 'siesta: Enegf =', DE_NEGF/eV,
352- . 'siesta: Exc =', Exc/eV,
353- . 'siesta: eta*DQ =', Ecorrec/eV,
354- . 'siesta: Emadel =', Emad/eV,
355- . 'siesta: Emeta =', Emeta/eV,
356- . 'siesta: Emolmec =', Emm/eV,
357- . 'siesta: Ekinion =', Ekinion/eV,
358- . 'siesta: Eharris =', (Eharrs+Ekinion)/eV,
359- . 'siesta: Etot =', (Etot+Ekinion)/eV,
360- . 'siesta: FreeEng =', (FreeE+Ekinion)/eV
361- if (cml_p) then
362- call cmlStartPropertyList(mainXML,
363- . title='Energy Decomposition')
364- call cmlAddProperty(xf=mainXML, value=Ebs/eV,
365- . units='siestaUnits:eV',
366- . dictref='siesta:Ebs', fmt='r6')
367- call cmlAddProperty(xf=mainXML, value=Eions/eV,
368- . units='siestaUnits:eV',
369- . dictref='siesta:Eions', fmt='r6')
370- call cmlAddProperty(xf=mainXML, value=Ena/eV,
371- . units='siestaUnits:eV',
372- . dictref='siesta:Ena', fmt='r6')
373- call cmlAddProperty(xf=mainXML, value=Ekin/eV,
374- . units='siestaUnits:eV',
375- . dictref='siesta:Ekin', fmt='r6')
376- call cmlAddProperty(xf=mainXML, value=Enl/eV,
377- . units='siestaUnits:eV',
378- . dictref='siesta:Enl', fmt='r6')
379- call cmlAddProperty(xf=mainXML, value=Eldau/eV,
380- . units='siestaUnits:eV',
381- . dictref='siesta:Eldau', fmt='r6')
382- call cmlAddProperty(xf=mainXML, value=DEna/eV,
383- . units='siestaUnits:eV',
384- . dictref='siesta:DEna', fmt='r6')
385- call cmlAddProperty(xf=mainXML, value=Eso/eV,
386- . units='siestaUnits:eV',
387- . dictref='siesta:Eso', fmt='r6')
388- call cmlAddProperty(xf=mainXML, value=DUscf/eV,
389- . units='siestaUnits:eV',
390- . dictref='siesta:DUscf', fmt='r6')
391- call cmlAddProperty(xf=mainXML, value=DUext/eV,
392- . units='siestaUnits:eV',
393- . dictref='siesta:DUext', fmt='r6')
394- call cmlAddProperty(xf=mainXML, value=DE_NEGF/eV,
395- . units='siestaUnits:eV',
396- . dictref='siesta:Enegf', fmt='r6')
397- call cmlAddProperty(xf=mainXML, value=Exc/eV,
398- . units='siestaUnits:eV',
399- . dictref='siesta:Exc', fmt='r6')
400- call cmlAddProperty(xf=mainXML,value=Ecorrec/eV,
401- . units='siestaUnits:eV',
402- . dictref='siesta:Ecorrec', fmt='r6')
403- call cmlAddProperty(xf=mainXML, value=Emad/eV,
404- . units='siestaUnits:eV',
405- . dictref='siesta:Emad', fmt='r6')
406- call cmlAddProperty(xf=mainXML, value=Emeta/eV,
407- . units='siestaUnits:eV',
408- . dictref='siesta:Emeta', fmt='r6')
409- call cmlAddProperty(xf=mainXML, value=Emm/eV,
410- . units='siestaUnits:eV',
411- . dictref='siesta:Emm', fmt='r6')
412- call cmlAddProperty(xf=mainXML,value=Ekinion/eV,
413- . units='siestaUnits:eV',
414- . dictref='siesta:Ekinion', fmt='r6')
415- call cmlAddProperty(xf=mainXML, value=(Eharrs+Ekinion)/eV,
416- . units='siestaUnits:eV',
417- . dictref='siesta:EharrsK', fmt='r6')
418- call cmlAddProperty(xf=mainXML, value=(Etot+Ekinion)/eV,
419- . units='siestaUnits:eV',
420- . dictref='siesta:EtotK', fmt='r6')
421- call cmlAddProperty(xf=mainXML, value=(FreeE+Ekinion)/eV,
422- . units='siestaUnits:eV',
423- . dictref='siesta:FreeEK', fmt='r6')
424- call cmlEndPropertyList(mainXML)
425- endif
426- endif
427-! On all SCF steps, print out the current energy
428- if (.not.final) then
429-! Print total energy and density matrix error
430- if (cml_p) then
431- call cmlStartPropertyList(mainXML, title='SCF Cycle')
432-! Eharrs is always output
433- call cmlAddProperty(xf=mainXML, value=Eharrs/eV,
434- . units="siestaUnits:eV",
435- . dictRef="siesta:Eharrs", fmt="r7")
436- if ( .not. harrisfun ) then
437- ! store dDmax and dHmax
438- call cmlAddProperty(xf=mainXML, value=dDmax,
439- . units="siestaUnits:none",
440- . dictRef="siesta:dDmax", fmt="r7")
441- call cmlAddProperty(xf=mainXML, value=dHmax/eV,
442- . units="siestaUnits:eV",
443- . dictRef="siesta:dHmax", fmt="r7")
444- end if
445- endif
446-
447-! Determines which properties are output.
448- if (harrisfun) then
449- write(6,"(/a,f14.6,/)") 'siesta: Eharris(eV) = ', Eharrs/eV
450-! No need for further cml output
451- elseif ((isolve==SOLVE_DIAGON) .or. (isolve==SOLVE_PEXSI)
452- . .or. (isolve==SOLVE_TRANSI)) then
453- if (cml_p) then
454- call cmlAddProperty(xf=mainXML, value=Etot/eV,
455- . units="siestaUnits:eV",
456- . dictRef="siesta:Etot", fmt="r7")
457- call cmlAddProperty(xf=mainXML, value=FreeE/eV,
458- . units="siestaUnits:eV",
459- . dictRef="siesta:FreeE", fmt="r7")
460- if ( fixspin ) then
461- call cmlAddProperty(xf=mainXML, value=Efs(1)/eV,
462- . units="siestaUnits:eV",
463- . dictRef="siesta:Ef_UP", fmt="r7")
464- call cmlAddProperty(xf=mainXML, value=Efs(2)/eV,
465- . units="siestaUnits:eV",
466- . dictRef="siesta:Ef_DN", fmt="r7")
467- call cmlAddProperty(xf=mainXML, value=Efs(:)/eV,
468- . units="siestaUnits:eV",
469- . dictRef="siesta:Efs")
470- else
471- call cmlAddProperty(xf=mainXML,value=Ef/eV,
472- . units="siestaUnits:eV",
473- . dictRef="siesta:Ef", fmt="r7")
474- end if
475- endif
476-
477+module write_subs_positions
478+
479+ implicit none
480+
481+ private
482+ public :: siesta_write_positions
483+
484+contains
485+
486+ subroutine siesta_write_positions(moved)
487+ use siesta_geom
488+ use atomlist, only: elem, iza
489+ use siesta_cml
490+ use units, only: Ang
491+ use zmatrix, only: lUseZmatrix, write_canonical_ucell_and_Zmatrix
492+ use m_iostruct, only: write_struct
493+ use siesta_options, only: idyn
494+
495+ ! Whether the structure has already changed in response to forces and stresses
496+ logical, intent(in) :: moved
497+
498+ ! Write out structural information in "crystallography" format,
499+ ! and in "canonical" (see Zmatrix module) zmatrix format.
500+
501+ call write_struct( ucell, na_u, isa, iza, xa, moved)
502+ if ( lUseZmatrix .and. (idyn == 0) ) then
503+ if ( moved ) then
504+ call write_canonical_ucell_and_Zmatrix(filename="NEXT_ITER.UCELL.ZMATRIX")
505+ else
506+ call write_canonical_ucell_and_Zmatrix(filename="OUT.UCELL.ZMATRIX")
507+ end if
508+ end if
509+
510+ ! Note that this information should be written only at the
511+ ! time of processing the current geometry (in state_init)
512+
513+ if ( (.not. moved) .and. cml_p ) then
514+ call cmlAddMolecule(xf=mainXML, natoms=na_u, elements=elem, &
515+ atomRefs=cisa, coords=xa/Ang)
516+ call cmlAddLattice(xf=mainXML, cell=ucell/Ang, &
517+ units='siestaUnits:Ang', dictref='siesta:ucell')
518+ end if
519+
520+ end subroutine siesta_write_positions
521+
522+end module write_subs_positions
523+
524+module write_subs_pressure
525+
526+ implicit none
527+
528+ private
529+ public :: siesta_write_stress_pressure
530+
531+contains
532+
533+ subroutine siesta_write_stress_pressure()
534+ use precision, only: dp
535+ use m_steps, only: final
536+ use siesta_geom, only: ucell, na_u, isa, xa
537+ use intrinsic_missing, only: VNORM
538+
539+ character(len=32) :: shape
540+ integer :: i, nbcell
541+ real(dp) :: bcell(3,3)
542+
543+ ! This collects in bcell the nbcell number of linearly combinations
544+ ! of the periodic vectors
545+ call shaper(ucell, na_u, isa, xa, shape, nbcell, bcell)
546+
547+ ! Before we continue we normalize the bcell vectors
548+ ! This ensures that scalar projections contains the *length*
549+ ! of the projection
550+ do i = 1 , nbcell
551+ bcell(:,i) = bcell(:,i) / VNORM(bcell(:,i))
552+ end do
553+
554+ if ( final ) then
555+ call write_stress_pressure_final(nbcell, bcell)
556+ else
557+ call write_stress_pressure_not_final(nbcell, bcell)
558+ end if
559+
560+ end subroutine siesta_write_stress_pressure
561+
562+ subroutine write_stress_pressure_not_final(nbcell, bcell)
563+ use precision, only: dp
564+ use siesta_cml
565+ use units
566+
567+ use siesta_options, only: RemoveIntraMolecularPressure, idyn, tp, varcel
568+ use siesta_geom, only: volume_of_some_cell
569+ use m_iostruct, only: write_struct
570+ use m_energies, only: FreeE
571+ use m_stress, only: stress, kin_stress, mstress, tstress
572+
573+ integer, intent(in) :: nbcell ! dimensionality of the simulation
574+ real(dp), intent(in) :: bcell(3,nbcell)
575+ integer :: ix
576+
577+ real(dp):: Psol ! Pressure of "solid"
578+ real(dp):: Press ! Pressure
579+ real(dp):: ps(3,3) ! Auxiliary array
580+
581+ ! Stress tensor and pressure:
582+
583+ ! Write Voigt components of total stress tensor
584+ ps = stress + kin_stress
585+ write(6,'(/,a,6f12.2)') 'Stress-tensor-Voigt (kbar):', &
586+ (ps(ix,ix)/kbar,ix=1,3), &
587+ ps(1,2)/kbar, ps(2,3)/kbar, ps(1,3)/kbar
588+
589+ Press = - ((ps(1,1) + ps(2,2) + ps(3,3))/3.0_dp)
590+ write(6,"(a,f14.4)") "(Free)E + p*V (eV/cell)", &
591+ (FreeE + Press*volume_of_some_cell)/eV
592+
593+ if ( RemoveIntraMolecularPressure ) then
594+ ps = mstress + kin_stress
595+ write(6,'(/,a,6f12.2)') 'Inter-Molecular-Stress-Voigt (kbar):', &
596+ (ps(ix,ix)/kbar,ix=1,3), &
597+ ps(1,2)/kbar, ps(2,3)/kbar, ps(1,3)/kbar
598+
599+ Press = - ((ps(1,1) + ps(2,2) + ps(3,3))/3.0_dp)
600+ write(6,"(a,f14.4)") "(Free)E + p_inter_molec * V (eV/cell)", &
601+ (FreeE + Press*volume_of_some_cell)/eV
602+ end if
603+
604+ ! This use of the volume is OK, as it is called from state_analysis,
605+ ! before possibly changing the cell.
606+ ! Write "target enthalpy" (E + pV, where p is the *target* pressure)
607+ write(6,"(a,f14.4)") "Target enthalpy (eV/cell)", &
608+ (FreeE + tp*volume_of_some_cell)/eV
609+
610+ ! Write stress to CML file always
611+ if ( cml_p ) then
612+ call cmlAddProperty(xf=mainXML, value=stress*Ang**3, &
613+ dictref='siesta:stress', title='Stress', &
614+ units='siestaUnits:evpa3')
615+ end if
616+
617+ ! Output depends on dynamics option
618+ select case ( idyn )
619+ case ( 0:5, 8, 9, 10 )
620+
621+ if (idyn==0 .and. (.not.varcel)) then
622+ continue
623+ else
624+ write(6,'(/,a,3(/,a,3f12.6))') 'siesta: Stress tensor (static) (eV/Ang**3):', &
625+ ' ',stress(:,1)*Ang**3/eV, &
626+ ' ',stress(:,2)*Ang**3/eV, &
627+ ' ',stress(:,3)*Ang**3/eV
628+
629+ select case ( nbcell )
630+ case ( 1 )
631+ call calc_stress_1D(nbcell, bcell, stress, ps)
632+ write(6,'(/,a,3(/,a,3f15.6))') &
633+ 'siesta: 1D stress tensor (static) (eV/Ang):', &
634+ ' ',ps(:,1)*Ang/eV, &
635+ ' ',ps(:,2)*Ang/eV, &
636+ ' ',ps(:,3)*Ang/eV
637+ case ( 2 )
638+ call calc_stress_2D(nbcell, bcell, stress, ps)
639+ write(6,'(/,a,3(/,a,3f15.6))') &
640+ 'siesta: 2D stress tensor (static) (eV/Ang**2):', &
641+ ' ',ps(:,1)*Ang**2/eV, &
642+ ' ',ps(:,2)*Ang**2/eV, &
643+ ' ',ps(:,3)*Ang**2/eV
644+ end select
645+
646+ Psol = - (stress(1,1) + stress(2,2) + stress(3,3)) / 3.0_dp
647+ write(6,'(/,a,f20.8,a)') 'siesta: Pressure (static):', Psol/kBar, ' kBar'
648+
649+ if ( cml_p ) then
650+ call cmlAddProperty(xf=mainXML, value=Psol, &
651+ dictref='siesta:psol', title='Pressure (Static)', &
652+ units='siestaUnits:kBar')
653+ end if
654+
655+ write(6,'(/,a,3(/,a,3f12.6))') 'siesta: Stress tensor (total) (eV/Ang**3):', &
656+ ' ',tstress(:,1)*Ang**3/eV, &
657+ ' ',tstress(:,2)*Ang**3/eV, &
658+ ' ',tstress(:,3)*Ang**3/eV
659+
660+ select case ( nbcell )
661+ case ( 1 )
662+ call calc_stress_1D(nbcell, bcell, tstress, ps)
663+ write(6,'(/,a,3(/,a,3f15.6))') &
664+ 'siesta: 1D stress tensor (total) (eV/Ang):', &
665+ ' ',ps(:,1)*Ang/eV, &
666+ ' ',ps(:,2)*Ang/eV, &
667+ ' ',ps(:,3)*Ang/eV
668+ case ( 2 )
669+ call calc_stress_2D(nbcell, bcell, tstress, ps)
670+ write(6,'(/,a,3(/,a,3f15.6))') &
671+ 'siesta: 2D stress tensor (total) (eV/Ang**2):', &
672+ ' ',ps(:,1)*Ang**2/eV, &
673+ ' ',ps(:,2)*Ang**2/eV, &
674+ ' ',ps(:,3)*Ang**2/eV
675+ end select
676+
677+ Psol = - (tstress(1,1) + tstress(2,2) + tstress(3,3)) / 3.0_dp
678+ write(6,'(/,a,f20.8,a)') 'siesta: Pressure (total):', Psol/kBar, ' kBar'
679+
680+ if ( cml_p ) then
681+ call cmlAddProperty(xf=mainXML, value=tstress*Ang**3, &
682+ dictref='siesta:tstress', title='Total Stress', &
683+ units='siestaUnits:evpa3')
684+ call cmlAddProperty(xf=mainXML, value=Psol, &
685+ dictref='siesta:tpsol', title='Pressure (Total)', &
686+ units='siestaUnits:kBar')
687+ end if
688+
689+ if ( RemoveIntraMolecularPressure ) then
690+
691+ write(6,'(/,a,3(/,a,3f12.6))') 'siesta: Stress tensor (nonmol) (eV/Ang**3):', &
692+ ' ',mstress(:,1)*Ang**3/eV, &
693+ ' ',mstress(:,2)*Ang**3/eV, &
694+ ' ',mstress(:,3)*Ang**3/eV
695+
696+ Psol = - (mstress(1,1) + mstress(2,2) + mstress(3,3)) / 3.0_dp
697+ write(6,'(/,a,f20.8,a)') 'siesta: Pressure (nonmol):', Psol/kBar, ' kBar'
698+
699+ if (cml_p) then
700+ call cmlAddProperty(xf=mainXML, value=mstress*Ang**3, &
701+ dictref='siesta:mstress', &
702+ title='Stress tensor (normal)', &
703+ units='siestaUnits:evpa3')
704+ call cmlAddProperty(xf=mainXML, value=Psol, &
705+ dictref='siesta:pmol', title='Pressure (Nonmol)', &
706+ units='siestaUnits:kBar')
707+ end if
708+
709+ ps = mstress + kin_stress
710+ write(6,'(/,a,3(/,a,3f12.6))') 'siesta: Stress tensor (nonmol+kin) (eV/Ang**3):', &
711+ ' ',ps(:,1)*Ang**3/eV, &
712+ ' ',ps(:,2)*Ang**3/eV, &
713+ ' ',ps(:,3)*Ang**3/eV
714+
715+ Psol = - (ps(1,1) + ps(2,2) + ps(3,3)) / 3.0_dp
716+ write(6,'(/,a,f20.8,a)') 'siesta: Pressure (nonmol+kin):', Psol/kBar, ' kBar'
717+
718+ if ( cml_p ) then
719+ call cmlAddProperty(xf=mainXML, value=ps*Ang**3, &
720+ dictref='siesta:tmstress', &
721+ title='Stress tensor (nonmol+kin)', &
722+ units='siestaUnits:evpa3')
723+ call cmlAddProperty(xf=mainXML, value=Psol, &
724+ dictref='siesta:tpmol', title='Pressure (Nonmol+Kin)', &
725+ units='siestaUnits:kBar')
726+ end if
727+ end if ! Remove intramolecular pressure
728+
729+ end if !varcel
730+
731+ end select !idyn
732+
733+ end subroutine write_stress_pressure_not_final
734+
735+ subroutine write_stress_pressure_final(nbcell, bcell)
736+ use precision, only: dp
737+ use siesta_cml
738+ use units
739+
740+ use siesta_geom, only: volume_of_some_cell
741+ use m_iostruct, only: write_struct
742+ use m_stress, only: stress, mstress, cstress
743+
744+ integer, intent(in) :: nbcell ! dimensionality of the simulation
745+ real(dp), intent(in) :: bcell(3,nbcell)
746+
747+ real(dp):: Pmol ! Molecular pressure (discounting Virial term)
748+ real(dp):: Psol ! Pressure of "solid"
749+ real(dp) :: ps(3,3) ! auxiliary array
750+
751+
752+ ! AG: Possible BUG
753+ ! The volume here refers to the "old" cell, and
754+ ! might be out of date if the cell has changed.
755+
756+ ! Print stress tensor unconditionally
757+ write(6,'(/,a,3(/,a,3f12.6))') 'siesta: Stress tensor (static) (eV/Ang**3):', &
758+ 'siesta: ',stress(:,1)*Ang**3/eV, &
759+ 'siesta: ',stress(:,2)*Ang**3/eV, &
760+ 'siesta: ',stress(:,3)*Ang**3/eV
761+ if ( cml_p ) then
762+ call cmlAddProperty(xf=mainXML, value=stress*Ang**3/eV, &
763+ dictref='siesta:stress', units='siestaUnits:eV_Ang__3')
764+ end if
765+
766+ select case ( nbcell )
767+ case ( 1 )
768+ call calc_stress_1D(nbcell, bcell, stress, ps)
769+ write(6,'(/,a,3(/,a,3f15.6))') &
770+ 'siesta: 1D stress tensor (static) (eV/Ang):', &
771+ 'siesta: ',ps(:,1)*Ang/eV, &
772+ 'siesta: ',ps(:,2)*Ang/eV, &
773+ 'siesta: ',ps(:,3)*Ang/eV
774+ case ( 2 )
775+ call calc_stress_2D(nbcell, bcell, stress, ps)
776+ write(6,'(/,a,3(/,a,3f15.6))') &
777+ 'siesta: 2D stress tensor (static) (eV/Ang**2):', &
778+ 'siesta: ',ps(:,1)*Ang**2/eV, &
779+ 'siesta: ',ps(:,2)*Ang**2/eV, &
780+ 'siesta: ',ps(:,3)*Ang**2/eV
781+ end select
782+
783+ ! Print constrained stress tensor if different from unconstrained
784+ if ( any(cstress /= stress ) ) then
785+ write(6,'(/,a,3(/,a,3f12.6))') 'siesta: Constrained stress tensor (static) (eV/Ang**3):', &
786+ 'siesta: ',cstress(:,1)*Ang**3/eV, &
787+ 'siesta: ',cstress(:,2)*Ang**3/eV, &
788+ 'siesta: ',cstress(:,3)*Ang**3/eV
789+ if ( cml_p ) then
790+ call cmlAddProperty(xf=mainXML, value=cstress*Ang**3/eV, &
791+ dictref='siesta:cstress', &
792+ units='siestaUnits:eV_Ang__3')
793+ end if
794+
795+ select case ( nbcell )
796+ case ( 1 )
797+ call calc_stress_1D(nbcell, bcell, cstress, ps)
798+ write(6,'(/,a,3(/,a,3f15.6))') &
799+ 'siesta: 1D constrained stress tensor (static) (eV/Ang):', &
800+ 'siesta: ',ps(:,1)*Ang/eV, &
801+ 'siesta: ',ps(:,2)*Ang/eV, &
802+ 'siesta: ',ps(:,3)*Ang/eV
803+ case ( 2 )
804+ call calc_stress_2D(nbcell, bcell, cstress, ps)
805+ write(6,'(/,a,3(/,a,3f15.6))') &
806+ 'siesta: 2D constrained stress tensor (static) (eV/Ang**2):', &
807+ 'siesta: ',ps(:,1)*Ang**2/eV, &
808+ 'siesta: ',ps(:,2)*Ang**2/eV, &
809+ 'siesta: ',ps(:,3)*Ang**2/eV
810+ end select
811+
812+ end if
813+
814+ ! Find pressure
815+ Psol = - (( stress(1,1) + stress(2,2) + stress(3,3) )/3.0_dp)
816+ Pmol = - (( mstress(1,1) + mstress(2,2) + mstress(3,3) )/3.0_dp)
817+
818+ write(6,'(/,a,f18.6,a)') 'siesta: Cell volume =', volume_of_some_cell/Ang**3, ' Ang**3'
819+ write(6,'(2(/,a),2a20,a,3(/,a,2f20.8,a))') &
820+ 'siesta: Pressure (static):', &
821+ 'siesta: ','Solid', 'Molecule', ' Units', &
822+ 'siesta: ', Psol, Pmol, ' Ry/Bohr**3', &
823+ 'siesta: ', Psol*Ang**3/eV, Pmol*Ang**3/eV, ' eV/Ang**3', &
824+ 'siesta: ', Psol/kBar, Pmol/kBar, ' kBar'
825+ if ( cml_p ) then
826+ call cmlStartPropertyList(mainXML, title='Final Pressure')
827+ call cmlAddProperty(xf=mainXML, &
828+ value=volume_of_some_cell/Ang**3, &
829+ title='cell volume', dictref='siesta:cellvol', &
830+ units='siestaUnits:Ang__3')
831+ call cmlAddProperty(xf=mainXML, value=Psol/kBar, &
832+ title='Pressure of Solid', dictref='siesta:pressSol', &
833+ units='siestaUnits:kbar')
834+ call cmlAddProperty(xf=mainXML, value=Pmol/kBar, &
835+ title='Pressure of Molecule', dictref='siesta:pressMol', &
836+ units='siestaUnits:kbar')
837+ call cmlEndPropertyList(mainXML)
838+ end if
839+
840+ end subroutine write_stress_pressure_final
841+
842+ !< Transfer the stress from a chain calculation from eV / Ang ** 3 -> eV / Ang
843+ subroutine calc_stress_1D(nbcell, bcell, stress, stress_1d)
844+ use precision, only: dp
845+ ! TODO for varcel it isn't necessarily volume_of_some_cell == VOLUME(ucell)
846+ use siesta_geom, only: ucell, volume_of_some_cell
847+ use intrinsic_missing, only: VNORM
848+
849+ integer, intent(in) :: nbcell
850+ real(dp), intent(in) :: bcell(3,nbcell)
851+ real(dp), intent(in) :: stress(3,3)
852+ real(dp), intent(out) :: stress_1d(3,3)
853+
854+ integer :: i_periodic
855+ real(dp) :: vlen
856+
857+ if ( nbcell /= 1 ) call die('calc_stress_1D: error in computation')
858+
859+ ! Locate the unit-cell vector with the largest component along the periodic direction
860+ i_periodic = maxloc( abs(matmul(bcell(:,1), ucell)), 1)
861+ vlen = abs( dot_product(bcell(:,1), ucell(:,i_periodic)) )
862+
863+ stress_1d(:,:) = 0._dp
864+ stress_1d(i_periodic,i_periodic) = stress(i_periodic,i_periodic)
865+
866+ stress_1d(:,:) = stress_1d(:, :) * volume_of_some_cell / vlen
867+
868+ end subroutine calc_stress_1D
869+
870+ !< Transfer the stress from a chain calculation from eV / Ang ** 3 -> eV / Ang ** 2
871+ subroutine calc_stress_2D(nbcell, bcell, stress, stress_2d)
872+ use precision, only: dp
873+ ! TODO for varcel it isn't necessarily volume_of_some_cell == VOLUME(ucell)
874+ use siesta_geom, only: ucell
875+ use intrinsic_missing, only: VNORM
876+
877+ integer, intent(in) :: nbcell
878+ real(dp), intent(in) :: bcell(3,nbcell)
879+ real(dp), intent(in) :: stress(3,3)
880+ real(dp), intent(out) :: stress_2d(3,3)
881+
882+ integer :: i_vacuum
883+ real(dp) :: vcross(3), vlen
884+
885+ if ( nbcell /= 2 ) call die('calc_stress_2D: error in computation')
886+
887+ ! Locate the unit-cell vector with the largest component along the non-periodic direction
888+ ! The non-periodic direction is the cross-product of the periodic directions
889+ call cross(bcell(:,1), bcell(:,2), vcross)
890+ vcross(:) = vcross(:) / VNORM(vcross)
891+ i_vacuum = maxloc( abs(matmul(vcross, ucell)), 1)
892+ vlen = abs( dot_product(vcross, ucell(:,i_vacuum)) )
893+
894+ stress_2d(:,:) = stress(:,:)
895+ stress_2d(:, i_vacuum) = 0._dp
896+ stress_2d(i_vacuum, :) = 0._dp
897+
898+ stress_2d(:,:) = stress_2d(:,:) * vlen
899+
900+ end subroutine calc_stress_2D
901+
902+end module write_subs_pressure
903+
904+module write_subs_energies
905+
906+ implicit none
907+
908+ private
909+ public :: siesta_write_energies
910+
911+contains
912+
913+ subroutine siesta_write_energies( iscf, dDmax, dHmax )
914+ use m_steps, only: final, istp
915+
916+ use siesta_options
917+ use siesta_cml
918+ use units
919+ use m_energies
920+ use m_spin
921+ use m_ts_global_vars, only: TSrun
922+ implicit none
923+
924+ integer :: iscf
925+ real(dp), intent(in) :: dDmax ! Max. change in dens. matrix elem.
926+ real(dp), intent(in) :: dHmax ! Max. change in H elements
927+
928+ character(len=6) :: scf_name
929+
930+ logical :: first_scf_step
931+
932+ scf_name = 'scf'
933+ if ( TSrun ) then
934+ scf_name = 'ts-scf'
935+ end if
936+
937+ first_scf_step = (iscf == 1)
938+ ! Only print out full decomposition at very beginning and end.
939+ if ((istp==1.and.first_scf_step).or.final) then
940+ write(6,'(/,a,/,(a,f17.6))') &
941+ 'siesta: Program''s energy decomposition (eV):', &
942+ 'siesta: Ebs =', Ebs/eV, &
943+ 'siesta: Eions =', Eions/eV, &
944+ 'siesta: Ena =', Ena/eV, &
945+ 'siesta: Ekin =', Ekin/eV, &
946+ 'siesta: Enl =', Enl/eV, &
947+ 'siesta: Eso =', Eso/eV, &
948+ 'siesta: Eldau =', Eldau/eV, &
949+ 'siesta: DEna =', DEna/eV, &
950+ 'siesta: DUscf =', DUscf/eV, &
951+ 'siesta: DUext =', DUext/eV, &
952+ 'siesta: Enegf =', DE_NEGF/eV, &
953+ 'siesta: Exc =', Exc/eV, &
954+ 'siesta: eta*DQ =', Ecorrec/eV, &
955+ 'siesta: Emadel =', Emad/eV, &
956+ 'siesta: Emeta =', Emeta/eV, &
957+ 'siesta: Emolmec =', Emm/eV, &
958+ 'siesta: Ekinion =', Ekinion/eV, &
959+ 'siesta: Eharris =', (Eharrs+Ekinion)/eV, &
960+ 'siesta: Etot =', (Etot+Ekinion)/eV, &
961+ 'siesta: FreeEng =', (FreeE+Ekinion)/eV
962+
963+ if ( cml_p ) then
964+ call cmlStartPropertyList(mainXML, &
965+ title='Energy Decomposition')
966+ call cmlAddProperty(xf=mainXML, value=Ebs/eV, &
967+ units='siestaUnits:eV', &
968+ dictref='siesta:Ebs', fmt='r6')
969+ call cmlAddProperty(xf=mainXML, value=Eions/eV, &
970+ units='siestaUnits:eV', &
971+ dictref='siesta:Eions', fmt='r6')
972+ call cmlAddProperty(xf=mainXML, value=Ena/eV, &
973+ units='siestaUnits:eV', &
974+ dictref='siesta:Ena', fmt='r6')
975+ call cmlAddProperty(xf=mainXML, value=Ekin/eV, &
976+ units='siestaUnits:eV', &
977+ dictref='siesta:Ekin', fmt='r6')
978+ call cmlAddProperty(xf=mainXML, value=Enl/eV, &
979+ units='siestaUnits:eV', &
980+ dictref='siesta:Enl', fmt='r6')
981+ call cmlAddProperty(xf=mainXML, value=Eldau/eV, &
982+ units='siestaUnits:eV', &
983+ dictref='siesta:Eldau', fmt='r6')
984+ call cmlAddProperty(xf=mainXML, value=DEna/eV, &
985+ units='siestaUnits:eV', &
986+ dictref='siesta:DEna', fmt='r6')
987+ call cmlAddProperty(xf=mainXML, value=Eso/eV, &
988+ units='siestaUnits:eV', &
989+ dictref='siesta:Eso', fmt='r6')
990+ call cmlAddProperty(xf=mainXML, value=DUscf/eV, &
991+ units='siestaUnits:eV', &
992+ dictref='siesta:DUscf', fmt='r6')
993+ call cmlAddProperty(xf=mainXML, value=DUext/eV, &
994+ units='siestaUnits:eV', &
995+ dictref='siesta:DUext', fmt='r6')
996+ call cmlAddProperty(xf=mainXML, value=DE_NEGF/eV, &
997+ units='siestaUnits:eV', &
998+ dictref='siesta:Enegf', fmt='r6')
999+ call cmlAddProperty(xf=mainXML, value=Exc/eV, &
1000+ units='siestaUnits:eV', &
1001+ dictref='siesta:Exc', fmt='r6')
1002+ call cmlAddProperty(xf=mainXML,value=Ecorrec/eV, &
1003+ units='siestaUnits:eV', &
1004+ dictref='siesta:Ecorrec', fmt='r6')
1005+ call cmlAddProperty(xf=mainXML, value=Emad/eV, &
1006+ units='siestaUnits:eV', &
1007+ dictref='siesta:Emad', fmt='r6')
1008+ call cmlAddProperty(xf=mainXML, value=Emeta/eV, &
1009+ units='siestaUnits:eV', &
1010+ dictref='siesta:Emeta', fmt='r6')
1011+ call cmlAddProperty(xf=mainXML, value=Emm/eV, &
1012+ units='siestaUnits:eV', &
1013+ dictref='siesta:Emm', fmt='r6')
1014+ call cmlAddProperty(xf=mainXML,value=Ekinion/eV, &
1015+ units='siestaUnits:eV', &
1016+ dictref='siesta:Ekinion', fmt='r6')
1017+ call cmlAddProperty(xf=mainXML, value=(Eharrs+Ekinion)/eV, &
1018+ units='siestaUnits:eV', &
1019+ dictref='siesta:EharrsK', fmt='r6')
1020+ call cmlAddProperty(xf=mainXML, value=(Etot+Ekinion)/eV, &
1021+ units='siestaUnits:eV', &
1022+ dictref='siesta:EtotK', fmt='r6')
1023+ call cmlAddProperty(xf=mainXML, value=(FreeE+Ekinion)/eV, &
1024+ units='siestaUnits:eV', &
1025+ dictref='siesta:FreeEK', fmt='r6')
1026+ call cmlEndPropertyList(mainXML)
1027+ end if
1028+ end if
1029+
1030+ ! On all SCF steps, print out the current energy
1031+ if ( .not.final ) then
1032+ ! Print total energy and density matrix error
1033+ if ( cml_p ) then
1034+ call cmlStartPropertyList(mainXML, title='SCF Cycle')
1035+ ! Eharrs is always output
1036+ call cmlAddProperty(xf=mainXML, value=Eharrs/eV, &
1037+ units="siestaUnits:eV", &
1038+ dictRef="siesta:Eharrs", fmt="r7")
1039+ if ( .not. harrisfun ) then
1040+ ! store dDmax and dHmax
1041+ call cmlAddProperty(xf=mainXML, value=dDmax, &
1042+ units="siestaUnits:none", &
1043+ dictRef="siesta:dDmax", fmt="r7")
1044+ call cmlAddProperty(xf=mainXML, value=dHmax/eV, &
1045+ units="siestaUnits:eV", &
1046+ dictRef="siesta:dHmax", fmt="r7")
1047+ end if
1048+ end if
1049+
1050+ ! Determines which properties are output.
1051+ if ( harrisfun ) then
1052+ write(6,"(/a,f14.6,/)") 'siesta: Eharris(eV) = ', Eharrs/eV
1053+ ! No need for further cml output
1054+ else
1055+ select case ( isolve )
1056+ case ( SOLVE_DIAGON, SOLVE_PEXSI, SOLVE_TRANSI )
1057+
1058+ if ( cml_p ) then
1059+ call cmlAddProperty(xf=mainXML, value=Etot/eV, &
1060+ units="siestaUnits:eV", &
1061+ dictRef="siesta:Etot", fmt="r7")
1062+ call cmlAddProperty(xf=mainXML, value=FreeE/eV, &
1063+ units="siestaUnits:eV", &
1064+ dictRef="siesta:FreeE", fmt="r7")
1065+ if ( fixspin ) then
1066+ call cmlAddProperty(xf=mainXML, value=Efs(1)/eV, &
1067+ units="siestaUnits:eV", &
1068+ dictRef="siesta:Ef_UP", fmt="r7")
1069+ call cmlAddProperty(xf=mainXML, value=Efs(2)/eV, &
1070+ units="siestaUnits:eV", &
1071+ dictRef="siesta:Ef_DN", fmt="r7")
1072+ call cmlAddProperty(xf=mainXML, value=Efs(:)/eV, &
1073+ units="siestaUnits:eV", &
1074+ dictRef="siesta:Efs")
1075+ else
1076+ call cmlAddProperty(xf=mainXML,value=Ef/eV, &
1077+ units="siestaUnits:eV", &
1078+ dictRef="siesta:Ef", fmt="r7")
1079+ end if
1080+ end if
1081+
1082 if ( fixspin ) then
1083- if ( (iscf == 1) .or. muldeb ) then
1084- write(6,'(/,a12,3a16,4a10)')
1085- . 'iscf', 'Eharris(eV)', 'E_KS(eV)', 'FreeEng(eV)',
1086- . 'dDmax', 'Ef_up', 'Ef_dn(eV)','dHmax(eV)'
1087- end if
1088- write(6,'(a8,i4,3f16.6,4f10.6)')
1089- . trim(scf_name)//': ',iscf, Eharrs/eV, Etot/eV, FreeE/eV,
1090- . dDmax, Efs(1:2)/eV, dHmax/eV
1091- else ! fixspin
1092- if ( (iscf == 1) .or. muldeb ) then
1093- write(6,'(/,a12,3a16,3a10)')
1094- . 'iscf', 'Eharris(eV)', 'E_KS(eV)', 'FreeEng(eV)',
1095- . 'dDmax','Ef(eV)','dHmax(eV)'
1096- end if
1097- write(6,'(a8,i4,3f16.6,3f10.6)')
1098- . trim(scf_name)//': ',iscf, Eharrs/eV, Etot/eV, FreeE/eV,
1099- . dDmax, Ef /eV, dHmax/eV
1100- end if ! fixspin
1101-
1102- elseif ((isolve==SOLVE_ORDERN) .or. (isolve==SOLVE_MINIM)) then
1103+ if ( (iscf == 1) .or. muldeb ) then
1104+ write(6,'(/,a12,3a16,4a10)') &
1105+ 'iscf', 'Eharris(eV)', 'E_KS(eV)', 'FreeEng(eV)', &
1106+ 'dDmax', 'Ef_up', 'Ef_dn(eV)','dHmax(eV)'
1107+ end if
1108+ write(6,'(a8,i4,3f16.6,4f10.6)') &
1109+ trim(scf_name)//': ',iscf, Eharrs/eV, Etot/eV, FreeE/eV, &
1110+ dDmax, Efs(1:2)/eV, dHmax/eV
1111+ else ! fixspin
1112+ if ( (iscf == 1) .or. muldeb ) then
1113+ write(6,'(/,a12,3a16,3a10)') &
1114+ 'iscf', 'Eharris(eV)', 'E_KS(eV)', 'FreeEng(eV)', &
1115+ 'dDmax','Ef(eV)','dHmax(eV)'
1116+ end if
1117+ write(6,'(a8,i4,3f16.6,3f10.6)') &
1118+ trim(scf_name)//': ',iscf, Eharrs/eV, Etot/eV, FreeE/eV, &
1119+ dDmax, Ef/eV, dHmax/eV
1120+ end if ! fixspin
1121+
1122+ case ( SOLVE_ORDERN, SOLVE_MINIM )
1123
1124 write(6,'(/,a15,i4)') 'siesta: iscf = ',iscf
1125- write(6,'(a14,f15.4,a13,f15.4,a10,f10.6/)')
1126- . 'Eharris(eV) = ',Eharrs/eV,
1127- . 'E_KS(eV) = ',Etot/eV,'dDmax = ',dDmax
1128- if (cml_p) then
1129- call cmlAddProperty(xf=mainXML, value=Etot/eV,
1130- . units="siestaUnits:eV",
1131- . dictRef="siesta:Etot", fmt="r7")
1132- endif
1133-
1134- endif !harrisfun/isolve
1135+ write(6,'(a14,f15.4,a13,f15.4,a10,f10.6/)') &
1136+ 'Eharris(eV) = ',Eharrs/eV, &
1137+ 'E_KS(eV) = ',Etot/eV,'dDmax = ',dDmax
1138+ if ( cml_p ) then
1139+ call cmlAddProperty(xf=mainXML, value=Etot/eV, &
1140+ units="siestaUnits:eV", &
1141+ dictRef="siesta:Etot", fmt="r7")
1142+ end if
1143
1144- if (cml_p) then
1145- call cmlEndPropertyList(mainXML)
1146- endif
1147+ end select
1148+ end if !harrisfun/isolve
1149+
1150+ if ( cml_p ) then
1151+ call cmlEndPropertyList(mainXML)
1152+ end if
1153
1154- else !final
1155-! Print out additional information in finalization.
1156- write(6,'(/,a)') 'siesta: Final energy (eV):'
1157- write(6,'(a,a15,f15.6)')
1158- . 'siesta: ', 'Band Struct. =', Ebs/eV,
1159- . 'siesta: ', 'Kinetic =', Ekin/eV,
1160- . 'siesta: ', 'Hartree =', Uscf/eV,
1161- . 'siesta: ', 'Eldau =', Eldau/eV,
1162- . 'siesta: ', 'Eso =', Eso/eV,
1163- . 'siesta: ', 'Ext. field =', DUext/eV,
1164- . 'siesta: ', 'Enegf =', DE_NEGF/eV,
1165- . 'siesta: ', 'Exch.-corr. =', Exc/eV,
1166- . 'siesta: ', 'Ion-electron =', (Enascf+Enl+DUscf-Uscf-Uatm)/eV,
1167- . 'siesta: ', 'Ion-ion =', (Ena+Uatm-Enaatm-Eions)/eV,
1168- . 'siesta: ', 'Ekinion =', Ekinion/eV,
1169- . 'siesta: ', 'Total =', (Etot+Ekinion)/eV
1170- if ( fixspin ) then
1171- write(6,'(a,a15,f15.6)')
1172- . 'siesta: ', 'Fermi_up =', Efs(1)/eV,
1173- . 'siesta: ', 'Fermi_dn =', Efs(2)/eV
1174- else
1175- write(6,'(a,a15,f15.6)')
1176- . 'siesta: ', 'Fermi =', Ef/eV
1177- end if
1178- if (cml_p) then
1179- call cmlStartPropertyList(xf=mainXML, title='Final Energy')
1180- call cmlAddProperty(xf=mainXML, value=Ebs/eV,
1181- . units='siestaUnits:eV',
1182- . dictref='siesta:Ebs', fmt='r7')
1183- call cmlAddProperty(xf=mainXML, value=Ekin/eV,
1184- . units='siestaUnits:eV',
1185- . dictref='siesta:Ekin', fmt='r7')
1186- call cmlAddProperty(xf=mainXML, value=Uscf/eV,
1187- . units='siestaUnits:eV',
1188- . dictref='siesta:Uscf', fmt='r7')
1189- call cmlAddProperty(xf=mainXML, value=Eldau/eV,
1190- . units='siestaUnits:eV',
1191- . dictref='siesta:Eldau', fmt='r7')
1192- call cmlAddProperty(xf=mainXML, value=Eso/eV,
1193- . units='siestaUnits:eV',
1194- . dictref='siesta:Eso', fmt='r7')
1195- call cmlAddProperty(xf=mainXML, value=DUext/eV,
1196- . units='siestaUnits:eV',
1197- . dictref='siesta:DUext', fmt='r7')
1198- call cmlAddProperty(xf=mainXML, value=DE_NEGF/eV,
1199- . units='siestaUnits:eV',
1200- . dictref='siesta:Enegf', fmt='r7')
1201- call cmlAddProperty(xf=mainXML, value=Exc/eV,
1202- . units='siestaUnits:eV',
1203- . dictref='siesta:Exc', fmt='r7')
1204- call cmlAddProperty(xf=mainXML,
1205- . value=(Enascf+Enl+DUscf-Uscf-Uatm)/eV,
1206- . units='siestaUnits:eV',
1207- . dictref='siesta:I-e', fmt='r7')
1208- call cmlAddProperty(xf=mainXML,
1209- . value=(Ena+Uatm-Enaatm-Eions)/eV,
1210- . units='siestaUnits:eV',
1211- . dictref='siesta:I-I', fmt='r7')
1212- call cmlAddProperty(xf=mainXML, value=Ekinion/eV,
1213- . units='siestaUnits:eV',
1214- . dictref='siesta:Ekinion', fmt='r7')
1215- call cmlAddProperty(xf=mainXML, value=(Etot+Ekinion)/eV,
1216- . units='siestaUnits:eV',
1217- . dictref='siesta:Etot', fmt='r7')
1218- call cmlEndPropertyList(mainXML)
1219- endif !cml_p
1220- endif !final
1221-
1222- end subroutine siesta_write_energies
1223-
1224- END MODULE write_subs_energies
1225-!
1226-! Main module
1227-!
1228- Module write_subs
1229-
1230- use write_subs_pressure, only : siesta_write_stress_pressure
1231- use write_subs_energies, only : siesta_write_energies
1232- use write_subs_positions, only : siesta_write_positions
1233-
1234- private
1235- public :: siesta_write_forces, siesta_write_stress_pressure,
1236- & siesta_write_energies, siesta_write_positions
1237- CONTAINS
1238-
1239- subroutine siesta_write_forces(istep)
1240- use parallel, only: IOnode
1241- USE siesta_options
1242- use siesta_geom
1243- use siesta_cml
1244- use units, only : Ang, eV
1245- use m_forces
1246- use m_steps, only: final, inicoor
1247- use alloc, only : re_alloc, de_alloc
1248- implicit none
1249- integer, intent(in) :: istep
1250- integer :: ia, ix
1251- real(dp) :: cfmax, fmax, fres
1252- real(dp) :: ftot(3), cftot(3)
1253+ else !final
1254+
1255+ ! Print out additional information in finalization.
1256+ write(6,'(/,a)') 'siesta: Final energy (eV):'
1257+ write(6,'(a,a15,f15.6)') &
1258+ 'siesta: ', 'Band Struct. =', Ebs/eV, &
1259+ 'siesta: ', 'Kinetic =', Ekin/eV, &
1260+ 'siesta: ', 'Hartree =', Uscf/eV, &
1261+ 'siesta: ', 'Eldau =', Eldau/eV, &
1262+ 'siesta: ', 'Eso =', Eso/eV, &
1263+ 'siesta: ', 'Ext. field =', DUext/eV, &
1264+ 'siesta: ', 'Enegf =', DE_NEGF/eV, &
1265+ 'siesta: ', 'Exch.-corr. =', Exc/eV, &
1266+ 'siesta: ', 'Ion-electron =', (Enascf+Enl+DUscf-Uscf-Uatm)/eV, &
1267+ 'siesta: ', 'Ion-ion =', (Ena+Uatm-Enaatm-Eions)/eV, &
1268+ 'siesta: ', 'Ekinion =', Ekinion/eV, &
1269+ 'siesta: ', 'Total =', (Etot+Ekinion)/eV
1270+ if ( fixspin ) then
1271+ write(6,'(a,a15,f15.6)') &
1272+ 'siesta: ', 'Fermi_up =', Efs(1)/eV, &
1273+ 'siesta: ', 'Fermi_dn =', Efs(2)/eV
1274+ else
1275+ write(6,'(a,a15,f15.6)') &
1276+ 'siesta: ', 'Fermi =', Ef/eV
1277+ end if
1278+ if ( cml_p ) then
1279+ call cmlStartPropertyList(xf=mainXML, title='Final Energy')
1280+ call cmlAddProperty(xf=mainXML, value=Ebs/eV, &
1281+ units='siestaUnits:eV', &
1282+ dictref='siesta:Ebs', fmt='r7')
1283+ call cmlAddProperty(xf=mainXML, value=Ekin/eV, &
1284+ units='siestaUnits:eV', &
1285+ dictref='siesta:Ekin', fmt='r7')
1286+ call cmlAddProperty(xf=mainXML, value=Uscf/eV, &
1287+ units='siestaUnits:eV', &
1288+ dictref='siesta:Uscf', fmt='r7')
1289+ call cmlAddProperty(xf=mainXML, value=Eldau/eV, &
1290+ units='siestaUnits:eV', &
1291+ dictref='siesta:Eldau', fmt='r7')
1292+ call cmlAddProperty(xf=mainXML, value=Eso/eV, &
1293+ units='siestaUnits:eV', &
1294+ dictref='siesta:Eso', fmt='r7')
1295+ call cmlAddProperty(xf=mainXML, value=DUext/eV, &
1296+ units='siestaUnits:eV', &
1297+ dictref='siesta:DUext', fmt='r7')
1298+ call cmlAddProperty(xf=mainXML, value=DE_NEGF/eV, &
1299+ units='siestaUnits:eV', &
1300+ dictref='siesta:Enegf', fmt='r7')
1301+ call cmlAddProperty(xf=mainXML, value=Exc/eV, &
1302+ units='siestaUnits:eV', &
1303+ dictref='siesta:Exc', fmt='r7')
1304+ call cmlAddProperty(xf=mainXML, value=(Enascf+Enl+DUscf-Uscf-Uatm)/eV, &
1305+ units='siestaUnits:eV', &
1306+ dictref='siesta:I-e', fmt='r7')
1307+ call cmlAddProperty(xf=mainXML, value=(Ena+Uatm-Enaatm-Eions)/eV, &
1308+ units='siestaUnits:eV', &
1309+ dictref='siesta:I-I', fmt='r7')
1310+ call cmlAddProperty(xf=mainXML, value=Ekinion/eV, &
1311+ units='siestaUnits:eV', &
1312+ dictref='siesta:Ekinion', fmt='r7')
1313+ call cmlAddProperty(xf=mainXML, value=(Etot+Ekinion)/eV, &
1314+ units='siestaUnits:eV', &
1315+ dictref='siesta:Etot', fmt='r7')
1316+ call cmlEndPropertyList(mainXML)
1317+ end if !cml_p
1318+ end if !final
1319+
1320+ end subroutine siesta_write_energies
1321+
1322+end module write_subs_energies
1323+
1324+
1325+module write_subs_forces
1326+
1327+ implicit none
1328+
1329+ private
1330+ public :: siesta_write_forces
1331+
1332+contains
1333+
1334+ subroutine siesta_write_forces(istep)
1335+ use parallel, only: IOnode
1336+ use siesta_cml
1337+ use units
1338+
1339+ use siesta_options, only: dx, idyn, ftol, writef
1340+ use siesta_geom, only: na_u
1341+ use m_forces
1342+ use m_steps, only: final, inicoor
1343+
1344+ integer, intent(in) :: istep
1345+ integer :: ia, ix
1346+ real(dp) :: cfmax, fmax, fres
1347+ real(dp) :: ftot(3), cftot(3)
1348
1349 #ifdef DEBUG
1350- call write_debug( ' PRE siesta_write_forces' )
1351+ call write_debug( ' PRE siesta_write_forces' )
1352 #endif
1353- fmax = maxval(abs(fa))
1354- ftot = sum(fa, dim=2)
1355- cfmax = maxval(abs(cfa))
1356- cftot = sum(cfa, dim=2)
1357- fres = sqrt( sum(fa**2) / (na_u*3.0_dp))
1358-
1359- ! Almost the same forces output whether during simulation
1360- ! or at the end. Unfortunately not quite, therefore slightly
1361- ! tortuous logic below. If we are content to change format
1362- ! of output file slightly, this can be simplified.
1363- if (.not.final) then
1364- ! print forces to xml every step.
1365- ! output forces to stdout depending on writef
1366- if (cml_p) then
1367- call cmlStartPropertyList(mainXML, title='Forces')
1368+
1369+ ! TODO perhaps it would be better to *only* count max forces
1370+ ! etc for the non-ghost atoms? In this way it reflects the
1371+ ! actual geometry.
1372+ ! Yes, ghost atoms has zero force (due to their relatively high
1373+ ! mass), but it matters in the residue?
1374+
1375+ fmax = 0._dp
1376+ ftot(:) = 0._dp
1377+ cfmax = 0._dp
1378+ cftot(:) = 0._dp
1379+ fres = 0._dp
1380+ ! We might as well calculate all in one loop
1381+ ! This will be much more efficient for large quantities of atoms.
1382+ do ia = 1, na_u
1383+ ! TODO check for ghost atom here?, possibly skip and fix counter (na_u for fres @ end)
1384+ do ix = 1, 3
1385+ ftot(ix) = ftot(ix) + fa(ix,ia)
1386+ fmax = max(fmax, abs(fa(ix,ia)))
1387+ cftot(ix) = cftot(ix) + cfa(ix,ia)
1388+ cfmax = max(cfmax, abs(cfa(ix,ia)))
1389+ fres = fres + fa(ix,ia) ** 2
1390+ end do
1391+ end do
1392+ ! Finally calculate the residue
1393+ fres = sqrt( fres / real(na_u * 3, dp) )
1394+
1395+ ! Almost the same forces output whether during simulation
1396+ ! or at the end. Unfortunately not quite, therefore slightly
1397+ ! tortuous logic below. If we are content to change format
1398+ ! of output file slightly, this can be simplified.
1399+ if ( .not. final ) then
1400+ ! print forces to xml every step.
1401+ ! output forces to stdout depending on writef
1402+ if ( cml_p ) then
1403+ call cmlStartPropertyList(mainXML, title='Forces')
1404+ call cmlAddComment(mainXML,"Output: matrix fa(1:3,1:na_u)")
1405+ call cmlAddProperty(xf=mainXML, value=fa*Ang/eV, &
1406+ dictref='siesta:forces', units='siestaUnits:evpa')
1407+ call cmlAddProperty(xf=mainXML, value=ftot*Ang/eV, &
1408+ dictref='siesta:ftot', units='siestaUnits:evpa')
1409+ call cmlAddProperty(xf=mainXML, value=fmax*Ang/eV, &
1410+ dictref='siesta:fmax', units='siestaUnits:evpa')
1411+ call cmlAddProperty(xf=mainXML, value=fres*Ang/eV, &
1412+ dictref='siesta:fres', units='siestaUnits:evpa')
1413+ call cmlAddProperty(xf=mainXML, value=cfmax*Ang/eV, &
1414+ dictref='siesta:cfmax', units='siestaUnits:evpa')
1415+ call cmlEndPropertyList(mainXML)
1416+ end if
1417+
1418+ write(6,'(/,a)') 'siesta: Atomic forces (eV/Ang):'
1419+ if ( writef ) then
1420+ write(6,'(i6,3f12.6)')(ia,(fa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
1421+ end if
1422+
1423+ ! Always write forces to .FA file
1424+ call iofa( na_u, fa , .false.)
1425+ if ( any(fa/=cfa) ) then
1426+ ! if any one constraint is not the same
1427+ ! as the real forces, the FAC file
1428+ ! will also be created
1429+ call iofa( na_u, cfa , .true. )
1430+ end if
1431+
1432+ write(6,'(40("-"),/,a6,3f12.6)') 'Tot',(ftot(ix)*Ang/eV,ix=1,3)
1433+ write(6,'(40("-"),/,a6, f12.6)') 'Max',fmax*Ang/eV
1434+ write(6,'(a6,f12.6,a)')'Res',fres*Ang/eV, &
1435+ ' sqrt( Sum f_i^2 / 3N )'
1436+ write(6,'(40("-"),/,a6, f12.6,a)') 'Max',cfmax*Ang/eV, &
1437+ ' constrained'
1438+
1439+ ! Write Force Constant matrix if FC calculation ...
1440+ select case ( idyn )
1441+ case ( 6 )
1442+ ! If the istep is the first step, then it
1443+ ! must be the first
1444+ call ofc(fa,dx,na_u,.false.,istep==inicoor)
1445+ if ( any(fa/=cfa) ) then
1446+ call ofc(cfa,dx,na_u,.true.,istep==inicoor)
1447+ end if
1448+ case ( 7 )
1449+ if ( IOnode ) write(*,*) 'phonon support deactivated'
1450+ end select
1451+
1452+ else !not final
1453+
1454+ ! In finalization, only print forces if sufficiently large.
1455+ if ( fmax > ftol ) then
1456+ write(6,'(/,a)') 'siesta: Atomic forces (eV/Ang):'
1457+ write(6,'(a,i6,3f12.6)') &
1458+ ('siesta: ', ia,(fa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
1459+ write(6,'(a,40("-"),/,a,a6,3f12.6)') &
1460+ 'siesta: ','siesta: ','Tot',(ftot(ix)*Ang/eV,ix=1,3)
1461+
1462+ if ( cml_p ) then
1463+ call cmlStartPropertyList(mainXML, title='Force Summary')
1464 call cmlAddComment(mainXML,"Output: matrix fa(1:3,1:na_u)")
1465- call cmlAddProperty(xf=mainXML, value=fa*Ang/eV,
1466- . dictref='siesta:forces', units='siestaUnits:evpa')
1467- call cmlAddProperty(xf=mainXML, value=ftot*Ang/eV,
1468- . dictref='siesta:ftot', units='siestaUnits:evpa')
1469- call cmlAddProperty(xf=mainXML, value=fmax*Ang/eV,
1470- . dictref='siesta:fmax', units='siestaUnits:evpa')
1471- call cmlAddProperty(xf=mainXML, value=fres*Ang/eV,
1472- . dictref='siesta:fres', units='siestaUnits:evpa')
1473- call cmlAddProperty(xf=mainXML, value=cfmax*Ang/eV,
1474- . dictref='siesta:cfmax', units='siestaUnits:evpa')
1475+ call cmlAddProperty(xf=mainXML, value=fa*Ang/eV, &
1476+ dictref='siesta:forces', units='siestaUnits:evpa')
1477+ call cmlAddProperty(xf=mainXML, value=ftot*Ang/eV, &
1478+ dictref='siesta:ftot', units='siestaUnits:evpa')
1479 call cmlEndPropertyList(mainXML)
1480- endif
1481-
1482- write(6,'(/,a)') 'siesta: Atomic forces (eV/Ang):'
1483- if (writef) then
1484- write(6,'(i6,3f12.6)')(ia,(fa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
1485- endif
1486-
1487- ! Always write forces to .FA file
1488- call iofa( na_u, fa , .false.)
1489- if ( any(fa/=cfa) ) then
1490- ! if any one constraint is not the same
1491- ! as the real forces, the FAC file
1492- ! will also be created
1493- call iofa( na_u, cfa , .true. )
1494 end if
1495-
1496- write(6,'(40("-"),/,a6,3f12.6)') 'Tot',(ftot(ix)*Ang/eV,ix=1,3)
1497- write(6,'(40("-"),/,a6, f12.6)') 'Max',fmax*Ang/eV
1498- write(6,'(a6,f12.6,a)')'Res',fres*Ang/eV,
1499- . ' sqrt( Sum f_i^2 / 3N )'
1500- write(6,'(40("-"),/,a6, f12.6,a)') 'Max',cfmax*Ang/eV,
1501- . ' constrained'
1502-
1503- ! Write Force Constant matrix if FC calculation ...
1504- select case (idyn)
1505- case(6)
1506- ! If the istep is the first step, then it
1507- ! must be the first
1508- call ofc(fa,dx,na_u,.false.,istep==inicoor)
1509- if ( any(fa/=cfa) ) then
1510- call ofc(cfa,dx,na_u,.true.,istep==inicoor)
1511- end if
1512- case(7)
1513-! call phonon_write_forces(fa,na_u,ucell,istep)
1514- if (IOnode) write(*,*) 'phonon support deactivated'
1515- end select
1516-
1517- else !not final
1518-! In finalization, only print forces if sufficiently large.
1519- if (fmax .gt. ftol) then
1520- write(6,'(/,a)') 'siesta: Atomic forces (eV/Ang):'
1521- write(6,'(a,i6,3f12.6)')
1522- . ('siesta: ', ia,(fa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
1523- write(6,'(a,40("-"),/,a,a6,3f12.6)')
1524- . 'siesta: ','siesta: ','Tot',(ftot(ix)*Ang/eV,ix=1,3)
1525- if (cml_p) then
1526- call cmlStartPropertyList(mainXML, title='Force Summary')
1527- call cmlAddComment(mainXML,"Output: matrix fa(1:3,1:na_u)")
1528- call cmlAddProperty(xf=mainXML, value=fa*Ang/eV,
1529- . dictref='siesta:forces', units='siestaUnits:evpa')
1530- call cmlAddProperty(xf=mainXML, value=ftot*Ang/eV,
1531- . dictref='siesta:ftot', units='siestaUnits:evpa')
1532+
1533+ end if
1534+
1535+ if ( any(cfa /= fa) ) then
1536+ if ( cfmax > ftol ) then
1537+ write(6,'(/,a)') 'siesta: Constrained forces (eV/Ang):'
1538+ write(6,'(a,i6,3f12.6)') &
1539+ ('siesta: ',ia,(cfa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
1540+ write(6,'(a,40("-"),/,a,a4,3f12.6)') &
1541+ 'siesta: ','siesta: ','Tot',(cftot(ix)*Ang/eV,ix=1,3)
1542+
1543+ if ( cml_p ) then
1544+ call cmlStartPropertyList(mainXML, &
1545+ title='Constrained Force Summary')
1546+ call cmlAddComment(mainXML, &
1547+ "Output: matrix cfa(1:3,1:na_u)")
1548+ call cmlAddProperty(xf=mainXML, value=cfa*Ang/eV, &
1549+ dictref='siesta:cforces', units='siestaUnits:evpa')
1550+ call cmlAddProperty(xf=mainXML, value=cftot*Ang/eV, &
1551+ dictref='siesta:cftot', units='siestaUnits:evpa')
1552 call cmlEndPropertyList(mainXML)
1553- endif !cml_p
1554- endif
1555- if (Any(cfa /= fa)) then
1556- if (cfmax .gt. ftol) then
1557- write(6,'(/,a)') 'siesta: Constrained forces (eV/Ang):'
1558- write(6,'(a,i6,3f12.6)')
1559- . ('siesta: ',ia,(cfa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
1560- write(6,'(a,40("-"),/,a,a4,3f12.6)')
1561- . 'siesta: ','siesta: ','Tot',(cftot(ix)*Ang/eV,ix=1,3)
1562- if (cml_p) then
1563- call cmlStartPropertyList(mainXML,
1564- . title='Constrained Force Summary')
1565- call cmlAddComment(mainXML,
1566- $ "Output: matrix cfa(1:3,1:na_u)")
1567- call cmlAddProperty(xf=mainXML, value=cfa*Ang/eV,
1568- . dictref='siesta:cforces', units='siestaUnits:evpa')
1569- call cmlAddProperty(xf=mainXML, value=cftot*Ang/eV,
1570- . dictref='siesta:cftot', units='siestaUnits:evpa')
1571- call cmlEndPropertyList(mainXML)
1572- endif !cml_p
1573- endif
1574- endif
1575- endif !final for forces
1576+ end if
1577+
1578+ end if
1579+ end if
1580+
1581+ end if !final for forces
1582 #ifdef DEBUG
1583- call write_debug( ' POS siesta_write_forces' )
1584+ call write_debug( ' POS siesta_write_forces' )
1585 #endif
1586-
1587- end subroutine siesta_write_forces
1588-
1589- End Module write_subs
1590+
1591+ end subroutine siesta_write_forces
1592+
1593+end module write_subs_forces
1594+
1595+! Main module
1596+module write_subs
1597+
1598+ use write_subs_pressure, only : siesta_write_stress_pressure
1599+ use write_subs_energies, only : siesta_write_energies
1600+ use write_subs_positions, only : siesta_write_positions
1601+ use write_subs_forces, only : siesta_write_forces
1602+
1603+ implicit none
1604+
1605+ public
1606+ public :: siesta_write_stress_pressure
1607+ public :: siesta_write_energies
1608+ public :: siesta_write_positions
1609+ public :: siesta_write_forces
1610+
1611+end module write_subs
1612
1613=== modified file 'Tests/Makefile'
1614--- Tests/Makefile 2019-01-30 14:06:15 +0000
1615+++ Tests/Makefile 2019-02-28 08:28:01 +0000
1616@@ -62,7 +62,7 @@
1617 si2x1h-dipole-gcs si2x1h-quench si64 si64_coop si_bandpoints \
1618 si_coop sic-slab si_fatbands sih sih_fire sih-mrrr sih_op_broyden \
1619 sinw sinw_2 si-optical si_pdos_gamma si_pdos_kgrid \
1620- var_cell wannier zmatrix
1621+ var_cell wannier zmatrix graphene_perp graphene_non_perp
1622
1623 # Tests that won't work with more than a few processors
1624 tests_mpi_max2 = h2_bessel h_chain h_chain2 \
1625
1626=== added file 'Tests/Reference/graphene_non_perp.out'
1627--- Tests/Reference/graphene_non_perp.out 1970-01-01 00:00:00 +0000
1628+++ Tests/Reference/graphene_non_perp.out 2019-02-28 08:28:01 +0000
1629@@ -0,0 +1,521 @@
1630+Siesta Version : siesta-4.1-1067--low-D-stress-4
1631+Architecture : x86_64-linux-gcc
1632+Compiler version: GNU Fortran (GCC) 8.3.0
1633+Compiler flags : mpifort -m64 -fPIC -O3 -ftree-vectorize -fexpensive-optimizations -funroll-loops -fprefetch-loop-arrays -fgraphite -fipa-icf -fipa-pure-const -march=native -fschedule-fusion -fselective-scheduling -fstack-arrays -fmax-stack-var-size=10000 -fipa-sra -fipa-cp -fno-second-underscore -Wrealloc-lhs -Wrealloc-lhs-all
1634+PP flags : -DSIESTA__FLOOK -DSIESTA__METIS -DSIESTA__FFTW -DTS_PVT_METIS -DSIESTA__MRRR -DMPI -DFC_HAVE_FLUSH -DFC_HAVE_ABORT -DCDF -DGRID_DP -DMPI_TIMING -DTRANSIESTA_TIMING -DTBTRANS_TIMING -DSIESTA__MUMPS -DNCDF_PARALLEL -DNCDF -DNCDF_4 -DSCALAPACK_DEBUG -D_DIAG_WORK -DSIESTA__ELPA
1635+Libraries : libncdf.a libfdict.a -lzmumps -lmumps_common -lpord -lparmetis -lmetis -lesmumps -lscotch -lscotcherr -L/opt/gnu/8.3.0/elpa/2018.11.001/lib -Wl,-rpath=/opt/gnu/8.3.0/elpa/2018.11.001/lib -lelpa -L/opt/gnu/8.3.0/scalapack/204/lib -Wl,-rpath=/opt/gnu/8.3.0/scalapack/204/lib -lscalapack -L/opt/gnu/8.3.0/openblas/0.3.5/lib -Wl,-rpath=/opt/gnu/8.3.0/openblas/0.3.5/lib -lopenblas -L/home/nicpa/codes/flook/obj -lflookall -ldl -lnetcdff -lnetcdf -lpnetcdf -lhdf5_hl -lhdf5 -lz -m64 -fPIC -O3 -ftree-vectorize -fexpensive-optimizations -funroll-loops -fprefetch-loop-arrays -fgraphite -fipa-icf -fipa-pure-const -march=native -fschedule-fusion -fselective-scheduling -fstack-arrays -fmax-stack-var-size=10000 -fipa-sra -fipa-cp -fno-second-underscore -Wrealloc-lhs -Wrealloc-lhs-all -L/opt/gnu/8.3.0/fftw/3.3.8/lib -Wl,-rpath=/opt/gnu/8.3.0/fftw/3.3.8/lib -lfftw3
1636+PARALLEL version
1637+NetCDF support
1638+NetCDF-4 support
1639+NetCDF-4 MPI-IO support
1640+METIS ordering support
1641+Lua support
1642+
1643+* Running on 2 nodes in parallel
1644+>> Start of run: 28-FEB-2019 9:23:47
1645+
1646+ ***********************
1647+ * WELCOME TO SIESTA *
1648+ ***********************
1649+
1650+reinit: Reading from ../graphene_non_perp.fdf
1651+
1652+reinit: -----------------------------------------------------------------------
1653+reinit: System Name:
1654+reinit: -----------------------------------------------------------------------
1655+reinit: System Label: graphene_non_perp
1656+reinit: -----------------------------------------------------------------------
1657+
1658+initatom: Reading input for the pseudopotentials and atomic orbitals ----------
1659+Species number: 1 Atomic number: 6 Label: C
1660+
1661+Ground state valence configuration: 2s02 2p02
1662+Reading pseudopotential information in formatted form from C.psf
1663+
1664+Valence configuration for pseudopotential generation:
1665+2s( 2.00) rc: 1.49
1666+2p( 2.00) rc: 1.52
1667+3d( 0.00) rc: 1.58
1668+For C, standard SIESTA heuristics set lmxkb to 3
1669+ (one more than the basis l, including polarization orbitals).
1670+Use PS.lmax or PS.KBprojectors blocks to override.
1671+C pseudopotential only contains V_ls up to l=2 -- lmxkb reset.
1672+
1673+<basis_specs>
1674+===============================================================================
1675+C Z= 6 Mass= 12.010 Charge= 0.17977+309
1676+Lmxo=1 Lmxkb= 2 BasisType=split Semic=F
1677+L=0 Nsemic=0 Cnfigmx=2
1678+ i=1 nzeta=2 polorb=0 (2s)
1679+ splnorm: 0.15000
1680+ vcte: 0.0000
1681+ rinn: 0.0000
1682+ qcoe: 0.0000
1683+ qyuk: 0.0000
1684+ qwid: 0.10000E-01
1685+ rcs: 0.0000 0.0000
1686+ lambdas: 1.0000 1.0000
1687+L=1 Nsemic=0 Cnfigmx=2
1688+ i=1 nzeta=2 polorb=1 (2p)
1689+ splnorm: 0.15000
1690+ vcte: 0.0000
1691+ rinn: 0.0000
1692+ qcoe: 0.0000
1693+ qyuk: 0.0000
1694+ qwid: 0.10000E-01
1695+ rcs: 0.0000 0.0000
1696+ lambdas: 1.0000 1.0000
1697+-------------------------------------------------------------------------------
1698+L=0 Nkbl=1 erefs: 0.17977+309
1699+L=1 Nkbl=1 erefs: 0.17977+309
1700+L=2 Nkbl=1 erefs: 0.17977+309
1701+===============================================================================
1702+</basis_specs>
1703+
1704+Filter: Requested cutoff: 100.000000 Ry
1705+atom: Called for C (Z = 6)
1706+
1707+read_vps: Pseudopotential generation method:
1708+read_vps: ATM3 Troullier-Martins
1709+Total valence charge: 4.00000
1710+
1711+xc_check: Exchange-correlation functional:
1712+xc_check: Ceperley-Alder
1713+V l=0 = -2*Zval/r beyond r= 1.4666
1714+V l=1 = -2*Zval/r beyond r= 1.5038
1715+V l=2 = -2*Zval/r beyond r= 1.5612
1716+All V_l potentials equal beyond r= 1.5612
1717+This should be close to max(r_c) in ps generation
1718+All pots = -2*Zval/r beyond r= 1.5612
1719+
1720+VLOCAL1: 99.0% of the norm of Vloc inside 17.809 Ry
1721+VLOCAL1: 99.9% of the norm of Vloc inside 40.586 Ry
1722+atom: Maximum radius for 4*pi*r*r*local-pseudopot. charge 1.88329
1723+atom: Maximum radius for r*vlocal+2*Zval: 1.62091
1724+GHOST: No ghost state for L = 0
1725+GHOST: No ghost state for L = 1
1726+GHOST: No ghost state for L = 2
1727+
1728+KBgen: Kleinman-Bylander projectors:
1729+ l= 0 rc= 1.747182 el= -1.001947 Ekb= 5.181700 kbcos= 0.300603
1730+ l= 1 rc= 1.747182 el= -0.398598 Ekb= -4.328763 kbcos= -0.367074
1731+ l= 2 rc= 1.955272 el= 0.002326 Ekb= -1.016175 kbcos= -0.009979
1732+
1733+KBgen: Total number of Kleinman-Bylander projectors: 9
1734+atom: -------------------------------------------------------------------------
1735+
1736+atom: SANKEY-TYPE ORBITALS:
1737+atom: Selected multiple-zeta basis: split
1738+
1739+SPLIT: Orbitals with angular momentum L= 0
1740+
1741+SPLIT: Basis orbitals for state 2s
1742+
1743+SPLIT: PAO cut-off radius determined from an
1744+SPLIT: energy shift= 0.003675 Ry
1745+
1746+ izeta = 1
1747+ lambda = 1.000000
1748+ rc = 5.120028
1749+ energy = -0.998768
1750+ kinetic = 0.827586
1751+ potential(screened) = -1.826353
1752+ potential(ionic) = -5.392266
1753+PAO Filter: Cutoff used: 49.00 Ry
1754+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
1755+PAO Filter: Nominal cutoff (before filtering): 26.80 Ry
1756+Filter: Number of eigenfunctions of the
1757+ filtering kernel (total, used)= 33 10
1758+
1759+ izeta = 2
1760+ rmatch = 3.518811
1761+ splitnorm = 0.150000
1762+ energy = -0.872553
1763+ kinetic = 1.267903
1764+ potential(screened) = -2.140455
1765+ potential(ionic) = -5.964330
1766+PAO Filter: Cutoff used: 49.00 Ry
1767+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
1768+PAO Filter: Nominal cutoff (before filtering): 31.56 Ry
1769+Filter: Number of eigenfunctions of the
1770+ filtering kernel (total, used)= 26 6
1771+
1772+SPLIT: Orbitals with angular momentum L= 1
1773+
1774+SPLIT: Basis orbitals for state 2p
1775+
1776+SPLIT: PAO cut-off radius determined from an
1777+SPLIT: energy shift= 0.003675 Ry
1778+
1779+ izeta = 1
1780+ lambda = 1.000000
1781+ rc = 6.253707
1782+ energy = -0.395147
1783+ kinetic = 2.375452
1784+ potential(screened) = -2.770599
1785+ potential(ionic) = -6.196142
1786+PAO Filter: Cutoff used: 49.00 Ry
1787+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
1788+PAO Filter: Nominal cutoff (before filtering): 41.54 Ry
1789+Filter: Number of eigenfunctions of the
1790+ filtering kernel (total, used)= 38 12
1791+
1792+ izeta = 2
1793+ rmatch = 3.745781
1794+ splitnorm = 0.150000
1795+ energy = -0.266202
1796+ kinetic = 3.550016
1797+ potential(screened) = -3.816218
1798+ potential(ionic) = -7.645193
1799+PAO Filter: Cutoff used: 49.00 Ry
1800+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
1801+PAO Filter: Nominal cutoff (before filtering): 44.66 Ry
1802+Filter: Number of eigenfunctions of the
1803+ filtering kernel (total, used)= 27 6
1804+
1805+POLgen: Perturbative polarization orbital with L= 2
1806+
1807+POLgen: Polarization orbital for state 2p
1808+
1809+ izeta = 1
1810+ rc = 6.253707
1811+ energy = 1.011661
1812+ kinetic = 2.138329
1813+ potential(screened) = -1.126669
1814+ potential(ionic) = -3.864107
1815+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
1816+PAO Filter: Nominal cutoff (before filtering): 55.94 Ry
1817+Filter: Number of eigenfunctions of the
1818+ filtering kernel (total, used)= 38 11
1819+atom: Total number of Sankey-type orbitals: 13
1820+
1821+atm_pop: Valence configuration (for local Pseudopot. screening):
1822+ 2s( 2.00)
1823+ 2p( 2.00)
1824+Vna: chval, zval: 4.00000 4.00000
1825+
1826+Vna: Cut-off radius for the neutral-atom potential: 6.253707
1827+VNA Filter: diagnostic kinetic energy tol: 0.003000 Ry
1828+VNA Filter: Nominal cutoff (before filtering): 21.67 Ry
1829+VNA Filter: Cutoff used: 100.00 Ry
1830+Filter: Number of eigenfunctions of the
1831+ filtering kernel (total, used)= 51 18
1832+
1833+atom: _________________________________________________________________________
1834+
1835+prinput: Basis input ----------------------------------------------------------
1836+
1837+PAO.BasisType split
1838+
1839+%block ChemicalSpeciesLabel
1840+ 1 6 C # Species index, atomic number, species label
1841+%endblock ChemicalSpeciesLabel
1842+
1843+%block PAO.Basis # Define Basis set
1844+C 2 # Species label, number of l-shells
1845+ n=2 0 2 # n, l, Nzeta
1846+ 5.120 3.519
1847+ 1.000 1.000
1848+ n=2 1 2 P 1 # n, l, Nzeta, Polarization, NzetaPol
1849+ 6.254 3.746
1850+ 1.000 1.000
1851+%endblock PAO.Basis
1852+
1853+prinput: ----------------------------------------------------------------------
1854+
1855+Dumping basis to NetCDF file C.ion.nc
1856+coor: Atomic-coordinates input format = Cartesian coordinates
1857+coor: (in Angstroms)
1858+
1859+siesta: Atomic coordinates (Bohr) and species
1860+siesta: 0.00000 0.00000 9.44863 1 1
1861+siesta: 2.68341 0.00000 9.44863 1 2
1862+
1863+siesta: System type = slab
1864+
1865+initatomlists: Number of atoms, orbitals, and projectors: 2 26 18
1866+
1867+siesta: ******************** Simulation parameters ****************************
1868+siesta:
1869+siesta: The following are some of the parameters of the simulation.
1870+siesta: A complete list of the parameters used, including default values,
1871+siesta: can be found in file out.fdf
1872+siesta:
1873+redata: Spin configuration = none
1874+redata: Number of spin components = 1
1875+redata: Time-Reversal Symmetry = T
1876+redata: Spin-spiral = F
1877+redata: Long output = F
1878+redata: Number of Atomic Species = 1
1879+redata: Charge density info will appear in .RHO file
1880+redata: Write Mulliken Pop. = NO
1881+redata: Matel table size (NRTAB) = 1024
1882+redata: Mesh Cutoff = 200.0000 Ry
1883+redata: Net charge of the system = 0.0000 |e|
1884+redata: Min. number of SCF Iter = 0
1885+redata: Max. number of SCF Iter = 1000
1886+redata: SCF convergence failure will abort job
1887+redata: SCF mix quantity = Hamiltonian
1888+redata: Mix DM or H after convergence = F
1889+redata: Recompute H after scf cycle = F
1890+redata: Mix DM in first SCF step = T
1891+redata: Write Pulay info on disk = F
1892+redata: New DM Mixing Weight = 0.2500
1893+redata: New DM Occupancy tolerance = 0.000000000001
1894+redata: No kicks to SCF
1895+redata: DM Mixing Weight for Kicks = 0.5000
1896+redata: Require Harris convergence for SCF = F
1897+redata: Harris energy tolerance for SCF = 0.000100 eV
1898+redata: Require DM convergence for SCF = T
1899+redata: DM tolerance for SCF = 0.000100
1900+redata: Require EDM convergence for SCF = F
1901+redata: EDM tolerance for SCF = 0.001000 eV
1902+redata: Require H convergence for SCF = T
1903+redata: Hamiltonian tolerance for SCF = 0.001000 eV
1904+redata: Require (free) Energy convergence for SCF = F
1905+redata: (free) Energy tolerance for SCF = 0.000100 eV
1906+redata: Using Saved Data (generic) = F
1907+redata: Use continuation files for DM = F
1908+redata: Neglect nonoverlap interactions = F
1909+redata: Method of Calculation = Diagonalization
1910+redata: Electronic Temperature = 299.9869 K
1911+redata: Fix the spin of the system = F
1912+redata: Dynamics option = Single-point calculation
1913+mix.SCF: Pulay mixing = Pulay
1914+mix.SCF: Variant = stable
1915+mix.SCF: History steps = 2
1916+mix.SCF: Linear mixing weight = 0.250000
1917+mix.SCF: Mixing weight = 0.250000
1918+mix.SCF: SVD condition = 0.1000E-07
1919+redata: Save all siesta data in one NC = F
1920+redata: ***********************************************************************
1921+
1922+%block SCF.Mixers
1923+ Pulay
1924+%endblock SCF.Mixers
1925+
1926+%block SCF.Mixer.Pulay
1927+ # Mixing method
1928+ method pulay
1929+ variant stable
1930+
1931+ # Mixing options
1932+ weight 0.2500
1933+ weight.linear 0.2500
1934+ history 2
1935+%endblock SCF.Mixer.Pulay
1936+
1937+DM_history_depth set to one: no extrapolation allowed by default for geometry relaxation
1938+Size of DM history Fstack: 1
1939+Total number of electrons: 8.000000
1940+Total ionic charge: 8.000000
1941+
1942+* ProcessorY, Blocksize: 1 14
1943+
1944+
1945+* Orbital distribution balance (max,min): 14 12
1946+
1947+ Kpoints in: 91 . Kpoints trimmed: 85
1948+
1949+siesta: k-grid: Number of k-points = 85
1950+siesta: k-grid: Cutoff (effective) = 15.987 Ang
1951+siesta: k-grid: Supercell and displacements
1952+siesta: k-grid: 13 0 0 0.000
1953+siesta: k-grid: 0 13 0 0.000
1954+siesta: k-grid: 0 0 1 0.000
1955+
1956+diag: Algorithm = D&C
1957+diag: Parallel over k = F
1958+diag: Use parallel 2D distribution = F
1959+diag: Parallel block-size = 14
1960+diag: Parallel distribution = 1 x 2
1961+diag: Used triangular part = Lower
1962+diag: Absolute tolerance = 0.100E-15
1963+diag: Orthogonalization factor = 0.100E-05
1964+diag: Memory factor = 1.0000
1965+
1966+superc: Internal auxiliary supercell: 9 x 9 x 1 = 81
1967+superc: Number of atoms, orbitals, and projectors: 162 2106 1458
1968+
1969+
1970+ts: **************************************************************
1971+ts: Save H and S matrices = F
1972+ts: Save DM and EDM matrices = F
1973+ts: Fix Hartree potential = F
1974+ts: Only save the overlap matrix S = F
1975+ts: **************************************************************
1976+
1977+************************ Begin: TS CHECKS AND WARNINGS ************************
1978+************************ End: TS CHECKS AND WARNINGS **************************
1979+
1980+
1981+ ====================================
1982+ Single-point calculation
1983+ ====================================
1984+
1985+superc: Internal auxiliary supercell: 9 x 9 x 1 = 81
1986+superc: Number of atoms, orbitals, and projectors: 162 2106 1458
1987+
1988+outcell: Unit cell vectors (Ang):
1989+ 2.130000 1.229756 0.000000
1990+ 2.130000 -1.229756 0.000000
1991+ 14.000000 14.000000 40.000000
1992+
1993+outcell: Cell vector modules (Ang) : 2.459512 2.459512 44.631827
1994+outcell: Cell angles (23,13,12) (deg): 83.4071 64.6281 60.0000
1995+outcell: Cell volume (Ang**3) : 209.5504
1996+<dSpData1D:S at geom step 0
1997+ <sparsity:sparsity for geom step 0
1998+ nrows_g=26 nrows=14 sparsity=17.6302 nnzs=11918, refcount: 7>
1999+ <dData1D:(new from dSpData1D) n=11918, refcount: 1>
2000+refcount: 1>
2001+new_DM -- step: 1
2002+Initializing Density Matrix...
2003+DM filled with atomic data:
2004+<dSpData2D:DM initialized from atoms
2005+ <sparsity:sparsity for geom step 0
2006+ nrows_g=26 nrows=14 sparsity=17.6302 nnzs=11918, refcount: 8>
2007+ <dData2D:DM n=11918 m=1, refcount: 1>
2008+refcount: 1>
2009+No. of atoms with KB's overlaping orbs in proc 0. Max # of overlaps: 25 265
2010+New grid distribution: 1
2011+ 1 1: 9 1: 9 1: 90
2012+ 2 1: 9 1: 9 91: 180
2013+
2014+InitMesh: MESH = 18 x 18 x 360 = 116640
2015+InitMesh: (bp) = 9 x 9 x 180 = 14580
2016+InitMesh: Mesh cutoff (required, used) = 200.000 200.612 Ry
2017+ExtMesh (bp) on 0 = 73 x 69 x 150 = 755550
2018+New grid distribution: 2
2019+ 1 1: 9 1: 9 1: 23
2020+ 2 1: 9 1: 9 24: 180
2021+New grid distribution: 3
2022+ 1 1: 9 1: 9 1: 23
2023+ 2 1: 9 1: 9 24: 180
2024+Setting up quadratic distribution...
2025+ExtMesh (bp) on 0 = 73 x 69 x 83 = 418071
2026+PhiOnMesh: Number of (b)points on node 0 = 1863
2027+PhiOnMesh: nlist on node 0 = 111880
2028+cdiag-debug: jobz=V, algo= 1, Node= 1, work= 1120, rwork= 1171, iwork= 200
2029+cdiag-debug: jobz=V, algo= 1, Node= 0, work= 1120, rwork= 1327, iwork= 200
2030+
2031+stepf: Fermi-Dirac step function
2032+
2033+siesta: Program's energy decomposition (eV):
2034+siesta: Ebs = -98.469477
2035+siesta: Eions = 498.569969
2036+siesta: Ena = 103.286514
2037+siesta: Ekin = 226.364281
2038+siesta: Enl = -32.545918
2039+siesta: Eso = 0.000000
2040+siesta: Eldau = 0.000000
2041+siesta: DEna = -14.724389
2042+siesta: DUscf = 1.474012
2043+siesta: DUext = 0.000000
2044+siesta: Enegf = 0.000000
2045+siesta: Exc = -95.509239
2046+siesta: eta*DQ = 0.000000
2047+siesta: Emadel = 0.000000
2048+siesta: Emeta = 0.000000
2049+siesta: Emolmec = 0.000000
2050+siesta: Ekinion = 0.000000
2051+siesta: Eharris = -305.657707
2052+siesta: Etot = -310.224707
2053+siesta: FreeEng = -310.224707
2054+
2055+ iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV) dHmax(eV)
2056+ scf: 1 -305.657707 -310.224707 -310.224707 1.827758 -7.441513 5.344366
2057+timer: Routine,Calls,Time,% = IterSCF 1 0.756 41.00
2058+ scf: 2 -310.499135 -310.365788 -310.365788 0.042575 -6.297738 3.458708
2059+ scf: 3 -310.561539 -310.479768 -310.479768 0.077215 -4.665977 0.123112
2060+ scf: 4 -310.479579 -310.479678 -310.479678 0.002052 -4.742725 0.046775
2061+ scf: 5 -310.479850 -310.479766 -310.479766 0.000771 -4.754931 0.033564
2062+ scf: 6 -310.479926 -310.479859 -310.479859 0.001916 -4.785919 0.000836
2063+ scf: 7 -310.479859 -310.479859 -310.479859 0.000025 -4.785865 0.000525
2064+
2065+SCF Convergence by DM+H criterion
2066+max |DM_out - DM_in| : 0.0000249642
2067+max |H_out - H_in| (eV) : 0.0005252509
2068+SCF cycle converged after 7 iterations
2069+
2070+Using DM_out to compute the final energy and forces
2071+No. of atoms with KB's overlaping orbs in proc 0. Max # of overlaps: 25 265
2072+
2073+siesta: E_KS(eV) = -310.4799
2074+
2075+siesta: E_KS - E_eggbox = -310.4799
2076+
2077+siesta: Atomic forces (eV/Ang):
2078+----------------------------------------
2079+ Tot 0.001315 -0.001258 -0.000178
2080+----------------------------------------
2081+ Max 0.001231
2082+ Res 0.000678 sqrt( Sum f_i^2 / 3N )
2083+----------------------------------------
2084+ Max 0.001231 constrained
2085+
2086+Stress-tensor-Voigt (kbar): -5.13 -5.09 0.01 0.03 0.01 -0.00
2087+(Free)E + p*V (eV/cell) -310.0348
2088+Target enthalpy (eV/cell) -310.4799
2089+
2090+siesta: Program's energy decomposition (eV):
2091+siesta: Ebs = -110.839772
2092+siesta: Eions = 498.569969
2093+siesta: Ena = 103.286514
2094+siesta: Ekin = 218.669880
2095+siesta: Enl = -30.520375
2096+siesta: Eso = 0.000000
2097+siesta: Eldau = 0.000000
2098+siesta: DEna = -10.040502
2099+siesta: DUscf = 0.867939
2100+siesta: DUext = 0.000000
2101+siesta: Enegf = 0.000000
2102+siesta: Exc = -94.173347
2103+siesta: eta*DQ = 0.000000
2104+siesta: Emadel = 0.000000
2105+siesta: Emeta = 0.000000
2106+siesta: Emolmec = 0.000000
2107+siesta: Ekinion = 0.000000
2108+siesta: Eharris = -310.479859
2109+siesta: Etot = -310.479859
2110+siesta: FreeEng = -310.479859
2111+
2112+siesta: Final energy (eV):
2113+siesta: Band Struct. = -110.839772
2114+siesta: Kinetic = 218.669880
2115+siesta: Hartree = 3382.532169
2116+siesta: Eldau = 0.000000
2117+siesta: Eso = 0.000000
2118+siesta: Ext. field = 0.000000
2119+siesta: Enegf = 0.000000
2120+siesta: Exch.-corr. = -94.173347
2121+siesta: Ion-electron = -6971.449101
2122+siesta: Ion-ion = 3153.940539
2123+siesta: Ekinion = 0.000000
2124+siesta: Total = -310.479859
2125+siesta: Fermi = -4.785865
2126+
2127+siesta: Stress tensor (static) (eV/Ang**3):
2128+siesta: -0.003199 0.000016 -0.000001
2129+siesta: 0.000016 -0.003177 0.000004
2130+siesta: -0.000001 0.000004 0.000006
2131+
2132+siesta: 2D stress tensor (static) (eV/Ang**2):
2133+siesta: -0.127979 0.000629 0.000000
2134+siesta: 0.000629 -0.127099 0.000000
2135+siesta: 0.000000 0.000000 0.000000
2136+
2137+siesta: Cell volume = 209.550434 Ang**3
2138+
2139+siesta: Pressure (static):
2140+siesta: Solid Molecule Units
2141+siesta: 0.00002313 0.00002314 Ry/Bohr**3
2142+siesta: 0.00212365 0.00212457 eV/Ang**3
2143+siesta: 3.40249403 3.40396886 kBar
2144+(Free)E+ p_basis*V_orbitals = -309.364883
2145+(Free)Eharris+ p_basis*V_orbitals = -309.364883
2146+
2147+siesta: Electric dipole (a.u.) = 0.000000 0.000000 -0.000003
2148+siesta: Electric dipole (Debye) = 0.000000 0.000000 -0.000008
2149+>> End of run: 28-FEB-2019 9:23:53
2150+Job completed
2151
2152=== added file 'Tests/Reference/graphene_perp.out'
2153--- Tests/Reference/graphene_perp.out 1970-01-01 00:00:00 +0000
2154+++ Tests/Reference/graphene_perp.out 2019-02-28 08:28:01 +0000
2155@@ -0,0 +1,521 @@
2156+Siesta Version : siesta-4.1-1067--low-D-stress-4
2157+Architecture : x86_64-linux-gcc
2158+Compiler version: GNU Fortran (GCC) 8.3.0
2159+Compiler flags : mpifort -m64 -fPIC -O3 -ftree-vectorize -fexpensive-optimizations -funroll-loops -fprefetch-loop-arrays -fgraphite -fipa-icf -fipa-pure-const -march=native -fschedule-fusion -fselective-scheduling -fstack-arrays -fmax-stack-var-size=10000 -fipa-sra -fipa-cp -fno-second-underscore -Wrealloc-lhs -Wrealloc-lhs-all
2160+PP flags : -DSIESTA__FLOOK -DSIESTA__METIS -DSIESTA__FFTW -DTS_PVT_METIS -DSIESTA__MRRR -DMPI -DFC_HAVE_FLUSH -DFC_HAVE_ABORT -DCDF -DGRID_DP -DMPI_TIMING -DTRANSIESTA_TIMING -DTBTRANS_TIMING -DSIESTA__MUMPS -DNCDF_PARALLEL -DNCDF -DNCDF_4 -DSCALAPACK_DEBUG -D_DIAG_WORK -DSIESTA__ELPA
2161+Libraries : libncdf.a libfdict.a -lzmumps -lmumps_common -lpord -lparmetis -lmetis -lesmumps -lscotch -lscotcherr -L/opt/gnu/8.3.0/elpa/2018.11.001/lib -Wl,-rpath=/opt/gnu/8.3.0/elpa/2018.11.001/lib -lelpa -L/opt/gnu/8.3.0/scalapack/204/lib -Wl,-rpath=/opt/gnu/8.3.0/scalapack/204/lib -lscalapack -L/opt/gnu/8.3.0/openblas/0.3.5/lib -Wl,-rpath=/opt/gnu/8.3.0/openblas/0.3.5/lib -lopenblas -L/home/nicpa/codes/flook/obj -lflookall -ldl -lnetcdff -lnetcdf -lpnetcdf -lhdf5_hl -lhdf5 -lz -m64 -fPIC -O3 -ftree-vectorize -fexpensive-optimizations -funroll-loops -fprefetch-loop-arrays -fgraphite -fipa-icf -fipa-pure-const -march=native -fschedule-fusion -fselective-scheduling -fstack-arrays -fmax-stack-var-size=10000 -fipa-sra -fipa-cp -fno-second-underscore -Wrealloc-lhs -Wrealloc-lhs-all -L/opt/gnu/8.3.0/fftw/3.3.8/lib -Wl,-rpath=/opt/gnu/8.3.0/fftw/3.3.8/lib -lfftw3
2162+PARALLEL version
2163+NetCDF support
2164+NetCDF-4 support
2165+NetCDF-4 MPI-IO support
2166+METIS ordering support
2167+Lua support
2168+
2169+* Running on 2 nodes in parallel
2170+>> Start of run: 28-FEB-2019 9:26:28
2171+
2172+ ***********************
2173+ * WELCOME TO SIESTA *
2174+ ***********************
2175+
2176+reinit: Reading from ../graphene_perp.fdf
2177+
2178+reinit: -----------------------------------------------------------------------
2179+reinit: System Name:
2180+reinit: -----------------------------------------------------------------------
2181+reinit: System Label: graphene_perp
2182+reinit: -----------------------------------------------------------------------
2183+
2184+initatom: Reading input for the pseudopotentials and atomic orbitals ----------
2185+Species number: 1 Atomic number: 6 Label: C
2186+
2187+Ground state valence configuration: 2s02 2p02
2188+Reading pseudopotential information in formatted form from C.psf
2189+
2190+Valence configuration for pseudopotential generation:
2191+2s( 2.00) rc: 1.49
2192+2p( 2.00) rc: 1.52
2193+3d( 0.00) rc: 1.58
2194+For C, standard SIESTA heuristics set lmxkb to 3
2195+ (one more than the basis l, including polarization orbitals).
2196+Use PS.lmax or PS.KBprojectors blocks to override.
2197+C pseudopotential only contains V_ls up to l=2 -- lmxkb reset.
2198+
2199+<basis_specs>
2200+===============================================================================
2201+C Z= 6 Mass= 12.010 Charge= 0.17977+309
2202+Lmxo=1 Lmxkb= 2 BasisType=split Semic=F
2203+L=0 Nsemic=0 Cnfigmx=2
2204+ i=1 nzeta=2 polorb=0 (2s)
2205+ splnorm: 0.15000
2206+ vcte: 0.0000
2207+ rinn: 0.0000
2208+ qcoe: 0.0000
2209+ qyuk: 0.0000
2210+ qwid: 0.10000E-01
2211+ rcs: 0.0000 0.0000
2212+ lambdas: 1.0000 1.0000
2213+L=1 Nsemic=0 Cnfigmx=2
2214+ i=1 nzeta=2 polorb=1 (2p)
2215+ splnorm: 0.15000
2216+ vcte: 0.0000
2217+ rinn: 0.0000
2218+ qcoe: 0.0000
2219+ qyuk: 0.0000
2220+ qwid: 0.10000E-01
2221+ rcs: 0.0000 0.0000
2222+ lambdas: 1.0000 1.0000
2223+-------------------------------------------------------------------------------
2224+L=0 Nkbl=1 erefs: 0.17977+309
2225+L=1 Nkbl=1 erefs: 0.17977+309
2226+L=2 Nkbl=1 erefs: 0.17977+309
2227+===============================================================================
2228+</basis_specs>
2229+
2230+Filter: Requested cutoff: 100.000000 Ry
2231+atom: Called for C (Z = 6)
2232+
2233+read_vps: Pseudopotential generation method:
2234+read_vps: ATM3 Troullier-Martins
2235+Total valence charge: 4.00000
2236+
2237+xc_check: Exchange-correlation functional:
2238+xc_check: Ceperley-Alder
2239+V l=0 = -2*Zval/r beyond r= 1.4666
2240+V l=1 = -2*Zval/r beyond r= 1.5038
2241+V l=2 = -2*Zval/r beyond r= 1.5612
2242+All V_l potentials equal beyond r= 1.5612
2243+This should be close to max(r_c) in ps generation
2244+All pots = -2*Zval/r beyond r= 1.5612
2245+
2246+VLOCAL1: 99.0% of the norm of Vloc inside 17.809 Ry
2247+VLOCAL1: 99.9% of the norm of Vloc inside 40.586 Ry
2248+atom: Maximum radius for 4*pi*r*r*local-pseudopot. charge 1.88329
2249+atom: Maximum radius for r*vlocal+2*Zval: 1.62091
2250+GHOST: No ghost state for L = 0
2251+GHOST: No ghost state for L = 1
2252+GHOST: No ghost state for L = 2
2253+
2254+KBgen: Kleinman-Bylander projectors:
2255+ l= 0 rc= 1.747182 el= -1.001947 Ekb= 5.181700 kbcos= 0.300603
2256+ l= 1 rc= 1.747182 el= -0.398598 Ekb= -4.328763 kbcos= -0.367074
2257+ l= 2 rc= 1.955272 el= 0.002326 Ekb= -1.016175 kbcos= -0.009979
2258+
2259+KBgen: Total number of Kleinman-Bylander projectors: 9
2260+atom: -------------------------------------------------------------------------
2261+
2262+atom: SANKEY-TYPE ORBITALS:
2263+atom: Selected multiple-zeta basis: split
2264+
2265+SPLIT: Orbitals with angular momentum L= 0
2266+
2267+SPLIT: Basis orbitals for state 2s
2268+
2269+SPLIT: PAO cut-off radius determined from an
2270+SPLIT: energy shift= 0.003675 Ry
2271+
2272+ izeta = 1
2273+ lambda = 1.000000
2274+ rc = 5.120028
2275+ energy = -0.998768
2276+ kinetic = 0.827586
2277+ potential(screened) = -1.826353
2278+ potential(ionic) = -5.392266
2279+PAO Filter: Cutoff used: 49.00 Ry
2280+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
2281+PAO Filter: Nominal cutoff (before filtering): 26.80 Ry
2282+Filter: Number of eigenfunctions of the
2283+ filtering kernel (total, used)= 33 10
2284+
2285+ izeta = 2
2286+ rmatch = 3.518811
2287+ splitnorm = 0.150000
2288+ energy = -0.872553
2289+ kinetic = 1.267903
2290+ potential(screened) = -2.140455
2291+ potential(ionic) = -5.964330
2292+PAO Filter: Cutoff used: 49.00 Ry
2293+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
2294+PAO Filter: Nominal cutoff (before filtering): 31.56 Ry
2295+Filter: Number of eigenfunctions of the
2296+ filtering kernel (total, used)= 26 6
2297+
2298+SPLIT: Orbitals with angular momentum L= 1
2299+
2300+SPLIT: Basis orbitals for state 2p
2301+
2302+SPLIT: PAO cut-off radius determined from an
2303+SPLIT: energy shift= 0.003675 Ry
2304+
2305+ izeta = 1
2306+ lambda = 1.000000
2307+ rc = 6.253707
2308+ energy = -0.395147
2309+ kinetic = 2.375452
2310+ potential(screened) = -2.770599
2311+ potential(ionic) = -6.196142
2312+PAO Filter: Cutoff used: 49.00 Ry
2313+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
2314+PAO Filter: Nominal cutoff (before filtering): 41.54 Ry
2315+Filter: Number of eigenfunctions of the
2316+ filtering kernel (total, used)= 38 12
2317+
2318+ izeta = 2
2319+ rmatch = 3.745781
2320+ splitnorm = 0.150000
2321+ energy = -0.266202
2322+ kinetic = 3.550016
2323+ potential(screened) = -3.816218
2324+ potential(ionic) = -7.645193
2325+PAO Filter: Cutoff used: 49.00 Ry
2326+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
2327+PAO Filter: Nominal cutoff (before filtering): 44.66 Ry
2328+Filter: Number of eigenfunctions of the
2329+ filtering kernel (total, used)= 27 6
2330+
2331+POLgen: Perturbative polarization orbital with L= 2
2332+
2333+POLgen: Polarization orbital for state 2p
2334+
2335+ izeta = 1
2336+ rc = 6.253707
2337+ energy = 1.011661
2338+ kinetic = 2.138329
2339+ potential(screened) = -1.126669
2340+ potential(ionic) = -3.864107
2341+PAO Filter: diagnostic kin energy tol: 0.003000 Ry
2342+PAO Filter: Nominal cutoff (before filtering): 55.94 Ry
2343+Filter: Number of eigenfunctions of the
2344+ filtering kernel (total, used)= 38 11
2345+atom: Total number of Sankey-type orbitals: 13
2346+
2347+atm_pop: Valence configuration (for local Pseudopot. screening):
2348+ 2s( 2.00)
2349+ 2p( 2.00)
2350+Vna: chval, zval: 4.00000 4.00000
2351+
2352+Vna: Cut-off radius for the neutral-atom potential: 6.253707
2353+VNA Filter: diagnostic kinetic energy tol: 0.003000 Ry
2354+VNA Filter: Nominal cutoff (before filtering): 21.67 Ry
2355+VNA Filter: Cutoff used: 100.00 Ry
2356+Filter: Number of eigenfunctions of the
2357+ filtering kernel (total, used)= 51 18
2358+
2359+atom: _________________________________________________________________________
2360+
2361+prinput: Basis input ----------------------------------------------------------
2362+
2363+PAO.BasisType split
2364+
2365+%block ChemicalSpeciesLabel
2366+ 1 6 C # Species index, atomic number, species label
2367+%endblock ChemicalSpeciesLabel
2368+
2369+%block PAO.Basis # Define Basis set
2370+C 2 # Species label, number of l-shells
2371+ n=2 0 2 # n, l, Nzeta
2372+ 5.120 3.519
2373+ 1.000 1.000
2374+ n=2 1 2 P 1 # n, l, Nzeta, Polarization, NzetaPol
2375+ 6.254 3.746
2376+ 1.000 1.000
2377+%endblock PAO.Basis
2378+
2379+prinput: ----------------------------------------------------------------------
2380+
2381+Dumping basis to NetCDF file C.ion.nc
2382+coor: Atomic-coordinates input format = Cartesian coordinates
2383+coor: (in Angstroms)
2384+
2385+siesta: Atomic coordinates (Bohr) and species
2386+siesta: 0.00000 0.00000 9.44863 1 1
2387+siesta: 2.68341 0.00000 9.44863 1 2
2388+
2389+siesta: System type = slab
2390+
2391+initatomlists: Number of atoms, orbitals, and projectors: 2 26 18
2392+
2393+siesta: ******************** Simulation parameters ****************************
2394+siesta:
2395+siesta: The following are some of the parameters of the simulation.
2396+siesta: A complete list of the parameters used, including default values,
2397+siesta: can be found in file out.fdf
2398+siesta:
2399+redata: Spin configuration = none
2400+redata: Number of spin components = 1
2401+redata: Time-Reversal Symmetry = T
2402+redata: Spin-spiral = F
2403+redata: Long output = F
2404+redata: Number of Atomic Species = 1
2405+redata: Charge density info will appear in .RHO file
2406+redata: Write Mulliken Pop. = NO
2407+redata: Matel table size (NRTAB) = 1024
2408+redata: Mesh Cutoff = 200.0000 Ry
2409+redata: Net charge of the system = 0.0000 |e|
2410+redata: Min. number of SCF Iter = 0
2411+redata: Max. number of SCF Iter = 1000
2412+redata: SCF convergence failure will abort job
2413+redata: SCF mix quantity = Hamiltonian
2414+redata: Mix DM or H after convergence = F
2415+redata: Recompute H after scf cycle = F
2416+redata: Mix DM in first SCF step = T
2417+redata: Write Pulay info on disk = F
2418+redata: New DM Mixing Weight = 0.2500
2419+redata: New DM Occupancy tolerance = 0.000000000001
2420+redata: No kicks to SCF
2421+redata: DM Mixing Weight for Kicks = 0.5000
2422+redata: Require Harris convergence for SCF = F
2423+redata: Harris energy tolerance for SCF = 0.000100 eV
2424+redata: Require DM convergence for SCF = T
2425+redata: DM tolerance for SCF = 0.000100
2426+redata: Require EDM convergence for SCF = F
2427+redata: EDM tolerance for SCF = 0.001000 eV
2428+redata: Require H convergence for SCF = T
2429+redata: Hamiltonian tolerance for SCF = 0.001000 eV
2430+redata: Require (free) Energy convergence for SCF = F
2431+redata: (free) Energy tolerance for SCF = 0.000100 eV
2432+redata: Using Saved Data (generic) = F
2433+redata: Use continuation files for DM = F
2434+redata: Neglect nonoverlap interactions = F
2435+redata: Method of Calculation = Diagonalization
2436+redata: Electronic Temperature = 299.9869 K
2437+redata: Fix the spin of the system = F
2438+redata: Dynamics option = Single-point calculation
2439+mix.SCF: Pulay mixing = Pulay
2440+mix.SCF: Variant = stable
2441+mix.SCF: History steps = 2
2442+mix.SCF: Linear mixing weight = 0.250000
2443+mix.SCF: Mixing weight = 0.250000
2444+mix.SCF: SVD condition = 0.1000E-07
2445+redata: Save all siesta data in one NC = F
2446+redata: ***********************************************************************
2447+
2448+%block SCF.Mixers
2449+ Pulay
2450+%endblock SCF.Mixers
2451+
2452+%block SCF.Mixer.Pulay
2453+ # Mixing method
2454+ method pulay
2455+ variant stable
2456+
2457+ # Mixing options
2458+ weight 0.2500
2459+ weight.linear 0.2500
2460+ history 2
2461+%endblock SCF.Mixer.Pulay
2462+
2463+DM_history_depth set to one: no extrapolation allowed by default for geometry relaxation
2464+Size of DM history Fstack: 1
2465+Total number of electrons: 8.000000
2466+Total ionic charge: 8.000000
2467+
2468+* ProcessorY, Blocksize: 1 14
2469+
2470+
2471+* Orbital distribution balance (max,min): 14 12
2472+
2473+ Kpoints in: 91 . Kpoints trimmed: 85
2474+
2475+siesta: k-grid: Number of k-points = 85
2476+siesta: k-grid: Cutoff (effective) = 15.987 Ang
2477+siesta: k-grid: Supercell and displacements
2478+siesta: k-grid: 13 0 0 0.000
2479+siesta: k-grid: 0 13 0 0.000
2480+siesta: k-grid: 0 0 1 0.000
2481+
2482+diag: Algorithm = D&C
2483+diag: Parallel over k = F
2484+diag: Use parallel 2D distribution = F
2485+diag: Parallel block-size = 14
2486+diag: Parallel distribution = 1 x 2
2487+diag: Used triangular part = Lower
2488+diag: Absolute tolerance = 0.100E-15
2489+diag: Orthogonalization factor = 0.100E-05
2490+diag: Memory factor = 1.0000
2491+
2492+superc: Internal auxiliary supercell: 9 x 9 x 1 = 81
2493+superc: Number of atoms, orbitals, and projectors: 162 2106 1458
2494+
2495+
2496+ts: **************************************************************
2497+ts: Save H and S matrices = F
2498+ts: Save DM and EDM matrices = F
2499+ts: Fix Hartree potential = F
2500+ts: Only save the overlap matrix S = F
2501+ts: **************************************************************
2502+
2503+************************ Begin: TS CHECKS AND WARNINGS ************************
2504+************************ End: TS CHECKS AND WARNINGS **************************
2505+
2506+
2507+ ====================================
2508+ Single-point calculation
2509+ ====================================
2510+
2511+superc: Internal auxiliary supercell: 9 x 9 x 1 = 81
2512+superc: Number of atoms, orbitals, and projectors: 162 2106 1458
2513+
2514+outcell: Unit cell vectors (Ang):
2515+ 2.130000 1.229756 0.000000
2516+ 2.130000 -1.229756 0.000000
2517+ 0.000000 0.000000 40.000000
2518+
2519+outcell: Cell vector modules (Ang) : 2.459512 2.459512 40.000000
2520+outcell: Cell angles (23,13,12) (deg): 90.0000 90.0000 60.0000
2521+outcell: Cell volume (Ang**3) : 209.5504
2522+<dSpData1D:S at geom step 0
2523+ <sparsity:sparsity for geom step 0
2524+ nrows_g=26 nrows=14 sparsity=17.6302 nnzs=11918, refcount: 7>
2525+ <dData1D:(new from dSpData1D) n=11918, refcount: 1>
2526+refcount: 1>
2527+new_DM -- step: 1
2528+Initializing Density Matrix...
2529+DM filled with atomic data:
2530+<dSpData2D:DM initialized from atoms
2531+ <sparsity:sparsity for geom step 0
2532+ nrows_g=26 nrows=14 sparsity=17.6302 nnzs=11918, refcount: 8>
2533+ <dData2D:DM n=11918 m=1, refcount: 1>
2534+refcount: 1>
2535+No. of atoms with KB's overlaping orbs in proc 0. Max # of overlaps: 25 265
2536+New grid distribution: 1
2537+ 1 1: 10 1: 10 1: 90
2538+ 2 1: 10 1: 10 91: 180
2539+
2540+InitMesh: MESH = 20 x 20 x 360 = 144000
2541+InitMesh: (bp) = 10 x 10 x 180 = 18000
2542+InitMesh: Mesh cutoff (required, used) = 200.000 223.865 Ry
2543+ExtMesh (bp) on 0 = 74 x 74 x 150 = 821400
2544+New grid distribution: 2
2545+ 1 1: 10 1: 10 1: 23
2546+ 2 1: 10 1: 10 24: 180
2547+New grid distribution: 3
2548+ 1 1: 10 1: 10 1: 23
2549+ 2 1: 10 1: 10 24: 180
2550+Setting up quadratic distribution...
2551+ExtMesh (bp) on 0 = 74 x 74 x 83 = 454508
2552+PhiOnMesh: Number of (b)points on node 0 = 2300
2553+PhiOnMesh: nlist on node 0 = 136545
2554+cdiag-debug: jobz=V, algo= 1, Node= 1, work= 1120, rwork= 1171, iwork= 200
2555+cdiag-debug: jobz=V, algo= 1, Node= 0, work= 1120, rwork= 1327, iwork= 200
2556+
2557+stepf: Fermi-Dirac step function
2558+
2559+siesta: Program's energy decomposition (eV):
2560+siesta: Ebs = -98.463721
2561+siesta: Eions = 498.569969
2562+siesta: Ena = 103.286514
2563+siesta: Ekin = 226.366689
2564+siesta: Enl = -32.546992
2565+siesta: Eso = 0.000000
2566+siesta: Eldau = 0.000000
2567+siesta: DEna = -14.725750
2568+siesta: DUscf = 1.474163
2569+siesta: DUext = 0.000000
2570+siesta: Enegf = 0.000000
2571+siesta: Exc = -95.509614
2572+siesta: eta*DQ = 0.000000
2573+siesta: Emadel = 0.000000
2574+siesta: Emeta = 0.000000
2575+siesta: Emolmec = 0.000000
2576+siesta: Ekinion = 0.000000
2577+siesta: Eharris = -305.658036
2578+siesta: Etot = -310.224959
2579+siesta: FreeEng = -310.224959
2580+
2581+ iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV) dHmax(eV)
2582+ scf: 1 -305.658036 -310.224959 -310.224959 1.827492 -7.441079 5.345135
2583+timer: Routine,Calls,Time,% = IterSCF 1 0.804 41.36
2584+ scf: 2 -310.499413 -310.366049 -310.366049 0.042706 -6.297139 3.459204
2585+ scf: 3 -310.561689 -310.479906 -310.479906 0.077462 -4.665061 0.123113
2586+ scf: 4 -310.479722 -310.479821 -310.479821 0.002052 -4.741810 0.046772
2587+ scf: 5 -310.479993 -310.479909 -310.479909 0.000773 -4.754016 0.033562
2588+ scf: 6 -310.480064 -310.479997 -310.479997 0.001919 -4.785004 0.000835
2589+ scf: 7 -310.479997 -310.479997 -310.479997 0.000025 -4.784949 0.000525
2590+
2591+SCF Convergence by DM+H criterion
2592+max |DM_out - DM_in| : 0.0000250158
2593+max |H_out - H_in| (eV) : 0.0005248258
2594+SCF cycle converged after 7 iterations
2595+
2596+Using DM_out to compute the final energy and forces
2597+No. of atoms with KB's overlaping orbs in proc 0. Max # of overlaps: 25 265
2598+
2599+siesta: E_KS(eV) = -310.4800
2600+
2601+siesta: E_KS - E_eggbox = -310.4800
2602+
2603+siesta: Atomic forces (eV/Ang):
2604+----------------------------------------
2605+ Tot -0.000000 -0.000000 -0.000000
2606+----------------------------------------
2607+ Max 0.000000
2608+ Res 0.000000 sqrt( Sum f_i^2 / 3N )
2609+----------------------------------------
2610+ Max 0.000000 constrained
2611+
2612+Stress-tensor-Voigt (kbar): -5.10 -5.10 0.01 0.00 -0.00 -0.00
2613+(Free)E + p*V (eV/cell) -310.0356
2614+Target enthalpy (eV/cell) -310.4800
2615+
2616+siesta: Program's energy decomposition (eV):
2617+siesta: Ebs = -110.835617
2618+siesta: Eions = 498.569969
2619+siesta: Ena = 103.286514
2620+siesta: Ekin = 218.671654
2621+siesta: Enl = -30.521252
2622+siesta: Eso = 0.000000
2623+siesta: Eldau = 0.000000
2624+siesta: DEna = -10.041489
2625+siesta: DUscf = 0.868009
2626+siesta: DUext = 0.000000
2627+siesta: Enegf = 0.000000
2628+siesta: Exc = -94.173465
2629+siesta: eta*DQ = 0.000000
2630+siesta: Emadel = 0.000000
2631+siesta: Emeta = 0.000000
2632+siesta: Emolmec = 0.000000
2633+siesta: Ekinion = 0.000000
2634+siesta: Eharris = -310.479997
2635+siesta: Etot = -310.479997
2636+siesta: FreeEng = -310.479997
2637+
2638+siesta: Final energy (eV):
2639+siesta: Band Struct. = -110.835617
2640+siesta: Kinetic = 218.671654
2641+siesta: Hartree = 3382.536543
2642+siesta: Eldau = 0.000000
2643+siesta: Eso = 0.000000
2644+siesta: Ext. field = 0.000000
2645+siesta: Enegf = 0.000000
2646+siesta: Exch.-corr. = -94.173465
2647+siesta: Ion-electron = -6971.455273
2648+siesta: Ion-ion = 3153.940544
2649+siesta: Ekinion = 0.000000
2650+siesta: Total = -310.479997
2651+siesta: Fermi = -4.784949
2652+
2653+siesta: Stress tensor (static) (eV/Ang**3):
2654+siesta: -0.003185 -0.000000 0.000000
2655+siesta: 0.000000 -0.003185 0.000000
2656+siesta: -0.000000 -0.000000 0.000008
2657+
2658+siesta: 2D stress tensor (static) (eV/Ang**2):
2659+siesta: -0.127413 -0.000000 0.000000
2660+siesta: 0.000000 -0.127414 0.000000
2661+siesta: 0.000000 0.000000 0.000000
2662+
2663+siesta: Cell volume = 209.550434 Ang**3
2664+
2665+siesta: Pressure (static):
2666+siesta: Solid Molecule Units
2667+siesta: 0.00002310 0.00002310 Ry/Bohr**3
2668+siesta: 0.00212089 0.00212089 eV/Ang**3
2669+siesta: 3.39808153 3.39808168 kBar
2670+(Free)E+ p_basis*V_orbitals = -309.365021
2671+(Free)Eharris+ p_basis*V_orbitals = -309.365021
2672+
2673+siesta: Electric dipole (a.u.) = 0.000000 0.000000 -0.000000
2674+siesta: Electric dipole (Debye) = 0.000000 0.000000 -0.000000
2675+>> End of run: 28-FEB-2019 9:26:34
2676+Job completed
2677
2678=== added directory 'Tests/graphene_non_perp'
2679=== added file 'Tests/graphene_non_perp/graphene_non_perp.fdf'
2680--- Tests/graphene_non_perp/graphene_non_perp.fdf 1970-01-01 00:00:00 +0000
2681+++ Tests/graphene_non_perp/graphene_non_perp.fdf 2019-02-28 08:28:01 +0000
2682@@ -0,0 +1,33 @@
2683+# Non perpendicular graphene cell
2684+
2685+SystemLabel graphene_non_perp
2686+
2687+NumberOfSpecies 1
2688+%block ChemicalSpeciesLabel
2689+ 1 6 C
2690+%endblock ChemicalSpeciesLabel
2691+
2692+LatticeConstant 1.0 Ang
2693+%block LatticeVectors
2694+ 2.13000000 1.22975607 0.00000000
2695+ 2.13000000 -1.22975607 0.00000000
2696+ 8.00000000 8.00000000 40.00000000
2697+%endblock LatticeVectors
2698+
2699+AtomicCoordinatesFormat Ang
2700+NumberOfAtoms 2
2701+%block AtomicCoordinatesAndAtomicSpecies
2702+ 0.00000000 0.00000000 5.00000000 1
2703+ 1.42000000 0.00000000 5.00000000 1
2704+%endblock AtomicCoordinatesAndAtomicSpecies
2705+
2706+PAO.EnergyShift 50. meV
2707+
2708+FilterCutoff 100. Ry
2709+MeshCutoff 200. Ry
2710+
2711+%block kgrid.MonkhorstPack
2712+ 13 0 0 0.
2713+ 0 13 0 0.
2714+ 0 0 1 0.
2715+%endblock
2716
2717=== added file 'Tests/graphene_non_perp/graphene_non_perp.pseudos'
2718--- Tests/graphene_non_perp/graphene_non_perp.pseudos 1970-01-01 00:00:00 +0000
2719+++ Tests/graphene_non_perp/graphene_non_perp.pseudos 2019-02-28 08:28:01 +0000
2720@@ -0,0 +1,1 @@
2721+C
2722
2723=== added file 'Tests/graphene_non_perp/makefile'
2724--- Tests/graphene_non_perp/makefile 1970-01-01 00:00:00 +0000
2725+++ Tests/graphene_non_perp/makefile 2019-02-28 08:28:01 +0000
2726@@ -0,0 +1,6 @@
2727+#
2728+# Single-test makefile
2729+#
2730+name=graphene_non_perp
2731+#
2732+include ../test.mk
2733
2734=== added directory 'Tests/graphene_perp'
2735=== added file 'Tests/graphene_perp/graphene_perp.fdf'
2736--- Tests/graphene_perp/graphene_perp.fdf 1970-01-01 00:00:00 +0000
2737+++ Tests/graphene_perp/graphene_perp.fdf 2019-02-28 08:28:01 +0000
2738@@ -0,0 +1,33 @@
2739+# Perpendicular graphene cell
2740+
2741+SystemLabel graphene_perp
2742+
2743+NumberOfSpecies 1
2744+%block ChemicalSpeciesLabel
2745+ 1 6 C
2746+%endblock ChemicalSpeciesLabel
2747+
2748+LatticeConstant 1.0 Ang
2749+%block LatticeVectors
2750+ 2.13000000 1.22975607 0.00000000
2751+ 2.13000000 -1.22975607 0.00000000
2752+ 0.00000000 0.00000000 40.00000000
2753+%endblock LatticeVectors
2754+
2755+AtomicCoordinatesFormat Ang
2756+NumberOfAtoms 2
2757+%block AtomicCoordinatesAndAtomicSpecies
2758+ 0.00000000 0.00000000 5.00000000 1
2759+ 1.42000000 0.00000000 5.00000000 1
2760+%endblock AtomicCoordinatesAndAtomicSpecies
2761+
2762+PAO.EnergyShift 50. meV
2763+
2764+FilterCutoff 100. Ry
2765+MeshCutoff 200. Ry
2766+
2767+%block kgrid.MonkhorstPack
2768+ 13 0 0 0.
2769+ 0 13 0 0.
2770+ 0 0 1 0.
2771+%endblock
2772
2773=== added file 'Tests/graphene_perp/graphene_perp.pseudos'
2774--- Tests/graphene_perp/graphene_perp.pseudos 1970-01-01 00:00:00 +0000
2775+++ Tests/graphene_perp/graphene_perp.pseudos 2019-02-28 08:28:01 +0000
2776@@ -0,0 +1,1 @@
2777+C
2778
2779=== added file 'Tests/graphene_perp/makefile'
2780--- Tests/graphene_perp/makefile 1970-01-01 00:00:00 +0000
2781+++ Tests/graphene_perp/makefile 2019-02-28 08:28:01 +0000
2782@@ -0,0 +1,6 @@
2783+#
2784+# Single-test makefile
2785+#
2786+name=graphene_perp
2787+#
2788+include ../test.mk
2789
2790=== modified file 'version.info'
2791--- version.info 2019-02-20 19:42:33 +0000
2792+++ version.info 2019-02-28 08:28:01 +0000
2793@@ -1,1 +1,5 @@
2794+<<<<<<< TREE
2795 siesta-4.1-1068
2796+=======
2797+siesta-4.1-1067--low-D-stress-5
2798+>>>>>>> MERGE-SOURCE

Subscribers

People subscribed via source and target branches