« back to all changes in this revision
Viewing changes to Src/siesta.F
- browse files at revision 228
- compare with another revision
- download diff
- download tarball
- view history from revision 228
- Committer: Alberto Garcia
- Date: 2007-11-21 11:45:49 UTC
- mfrom: (192.1.67)
- Revision ID: Arch-1:siesta@uam.es--2006%siesta-devel--reference--2.3--patch-1
Direct merge into master branch of initial BSC changes
The changes along the BSC branches, up to the end of the restructuring
of siesta.F and associated changes, have been merged into a direct
descendant of the main development line. The BSC work originally
started as a branch of siesta-devel--reference--2.1--patch-29. Development
along 2.1 continued, and a new continuation branch 2.3 has been created
specifically for this merge.
Main patches applied:
* ref@bsc--2007/siesta-bsc--master--2.1--base-0
tag of siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-7
* ref@bsc--2007/siesta-bsc--master--2.1--patch-1
Use of xalast in analysis routines. Exit of geometry loop
* ref@bsc--2007/siesta-bsc--master--2.1--patch-2
Creation of a module to hold the siesta options
* ref@bsc--2007/siesta-bsc--master--2.1--patch-3
New geometry module
* ref@bsc--2007/siesta-bsc--master--2.1--patch-4
New stub module for sparse matrices
* ref@bsc--2007/siesta-bsc--master--2.1--patch-5
More options for running tests
* ref@bsc--2007/siesta-bsc--master--2.1--patch-6
Encapsulation of k-point handling
* ref@bsc--2007/siesta-bsc--master--2.1--patch-7
Initialize iza in struct_init
* ref@bsc--2007/siesta-bsc--master--2.1--patch-8
Siesta_todo slimming by M. Quero
* ref@bsc--2007/siesta-bsc--master--2.1--patch-9
Fix import of no_l in born_charge
* ref@bsc--2007/siesta-bsc--master--2.1--patch-10
Some conversions to Fortran90 by M. Quero
* ref@bsc--2007/siesta-bsc--master--2.1--patch-11
Fourth session at the BSC
* ref@bsc--2007/siesta-bsc--master--2.1--patch-12
Clarification of the scope of the stress variables
* ref@bsc--2007/siesta-bsc--master--2.1--patch-13
Creation of siesta_forces
* ref@bsc--2007/siesta-bsc--master--2.1--patch-14
Replacement of some allocatables by pointers and automatics
* ref@bsc--2007/siesta-bsc--master--2.1--patch-15
New m_energies and m_steps modules. Back to old k-point behavior
* ref@bsc--2007/siesta-bsc--master--2.1--patch-16
Merge of removal of integer and real variables from siesta_todo
* ref@bsc--2007/siesta-bsc--master--2.1--patch-17
Fixes for troublesome bugs in reference code
* ref@bsc--2007/siesta-bsc--master--2.1--patch-18
Final cleanup of siesta_todo
* ref@bsc--2007/siesta-bsc--master--2.1--patch-19
Add character(len=*) routines to alloc.F90
* ref@bsc--2007/siesta-bsc--master--2.1--patch-20
Pointers in fixed and setspatial. si2x1h test added to bsc-Makefile
* ref@bsc--2007/siesta-bsc--master--2.1--patch-21
Re-organization of pulay module
* ref@bsc--2007/siesta-bsc--master--2.1--patch-22
Reorganization of hsparse/xijorb calls with new neighbor module
* ref@bsc--2007/siesta-bsc--master--2.1--patch-23
VPATH-aware compilation for multiple executable versions
* ref@bsc--2007/siesta-bsc--master--2.1--patch-24
Allocatables to pointers I
* ref@bsc--2007/siesta-bsc--master--2.1--patch-25
Allocatables to pointers II
* ref@bsc--2007/siesta-bsc--master--2.1--patch-26
Allocatables to pointers III
* ref@bsc--2007/siesta-bsc--master--2.1--patch-27
Allocatables to pointers IV -- new neighbor code + nspecies fix
* ref@bsc--2007/siesta-bsc--master--2.1--patch-28
Explicit array-ness in calls in initatom and cellxc
* ref@bsc--2007/siesta-bsc--master--2.1--patch-29
Fix wrong allocations in cellxc.F
* ref@bsc--2007/siesta-bsc--master--2.1--patch-30
Explicit array extents in initatom.f
* ref@bsc--2007/siesta-bsc--master--2.1--patch-31
Avoid shrinking of density-matrix arrays for extrapol.
* ref@bsc--2007/siesta-bsc--master--2.1--patch-32
Fix typo in state_init.F
* ref@bsc--2007/siesta-bsc--master--2.1--patch-33
Execute SCF loop when nscf=1
* ref@bsc--2007/siesta-bsc--master--2.1--patch-34
Clarify bounds of SCF loop in siesta_forces.
* ref@bsc--2007/siesta-bsc--master--2.2--base-0
tag of ref@bsc--2007/siesta-bsc--master--2.1--patch-31
* ref@bsc--2007/siesta-bsc--master--2.3--base-0
tag of ref@bsc--2007/siesta-bsc--master--2.2--base-0
* ref@bsc--2007/siesta-bsc--master--2.3--patch-1
Prepare CHANGES file
* ref@bsc--2007/siesta-bsc--master--2.3--patch-2
Merge patch-log for patch-30 from 2.1
* ref@bsc--2007/siesta-bsc--master--2.3--patch-3
New treatment of fractional atoms in VCA. Bug fix in lmxo
* ref@bsc--2007/siesta-bsc--master--2.3--patch-4
Re-enabling of kgrid update in variable-cell calculations
* ref@bsc--2007/siesta-bsc--master--2.3--patch-5
Fix typo in state_init.F
* ref@bsc--2007/siesta-bsc--master--2.3--patch-6
Merge filtering package by Eduardo Anglada
* ref@bsc--2007/siesta-bsc--master--2.3--patch-7
Add graphite_c6_full test for more realistic vdW test
* ref@bsc--2007/siesta-bsc--master--2.3--patch-8
Add patchlog for k-point fix
* ref@bsc--2007/siesta-bsc--master--2.3--patch-9
Fix of zmatrix code to deal with degenerate case
* ref@bsc--2007/siesta-bsc--master--2.3--patch-10
Units conversion in Util/Optical/optical.f. Scripts. Cosmetics
* ref@bsc--2007/siesta-bsc--master--2.3--patch-11
Implementation of basis_enthalpy calculation. Zmatrix dependency
* ref@bsc--2007/siesta-bsc--master--2.3--patch-12
Option to use fractional rc's for multiple zeta
* ref@bsc--2007/siesta-bsc--master--2.3--patch-13
Sync vpath changes. Update siesta.tex
* ref@bsc--2007/siesta-bsc--master--2.3--patch-14
Fix tag in MPI send/receive in mulliken
* ref@bsc--2007/siesta-bsc--master--2.3--patch-15
Option to use fractional rc's for multiple zeta
* ref@bsc--2007/siesta-bsc--master--2.3--patch-16
Merge of XML tester by Eduardo Anglada. Portability fixes
* ref@bsc--2007/siesta-bsc--master--2.3--patch-17
Fix marenostrum-mpi.make. Syntax in compare_m.f90
* ref@bsc--2007/siesta-bsc--master--2.3--patch-18
MareNostrum fixes. zdrot to blas. obj_setup. compare_m syntax
* ref@bsc--2007/siesta-bsc--master--2.3--patch-19
Sync to bsc--2.1. Change banner in CHANGES to bsc--2.3
* ref@bsc--2007/siesta-bsc--master--2.3--patch-20
Execute SCF loop when nscf=1
* ref@bsc--2007/siesta-bsc--master--2.3--patch-21
Refinements of XML tester. New ioncat program
* ref@bsc--2007/siesta-bsc--master--2.3--patch-22
Undef var in kgridinit, iohs MPI write, pdosg array bound, Origin shift
* ref@bsc--2007/siesta-bsc--master--2.3--patch-23
Zmatrix optimization enhancements. Sign change in MM stress.
* ref@bsc--2007/siesta-bsc--master--2.3--patch-24
Clarify bounds of SCF loop in siesta_forces.
* ref@bsc--2007/siesta-bsc--master--2.3--patch-25
Avoid MP-grid permutations in trivial gamma case
* ref@bsc--2007/siesta-bsc--master--2.3--patch-26
Sync to reference--2.1
* siesta@uam.es--2006/siesta-bsc--reference--2.1--base-0
tag of siesta@uam.es--2006/siesta-devel--reference--2.1--patch-29
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-1
First stage of siesta.F splitting at BSC
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-2
Created struct_init for initial geometry setup. New test force_2
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-3
Consolidate geometry updates at the end of loop
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-4
Initialize vol2 correctly in m_check_supercell.f
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-5
Kgrid setup streamlined. Bands. Proximity check. Hsparse allocation
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-6
Work by Manuel Quero before the meeting on Feb 7th
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-7
Moved the Born-effective-charge code
The changes along the BSC branches, up to the end of the restructuring
of siesta.F and associated changes, have been merged into a direct
descendant of the main development line. The BSC work originally
started as a branch of siesta-devel--reference--2.1--patch-29. Development
along 2.1 continued, and a new continuation branch 2.3 has been created
specifically for this merge.
Main patches applied:
* ref@bsc--2007/siesta-bsc--master--2.1--base-0
tag of siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-7
* ref@bsc--2007/siesta-bsc--master--2.1--patch-1
Use of xalast in analysis routines. Exit of geometry loop
* ref@bsc--2007/siesta-bsc--master--2.1--patch-2
Creation of a module to hold the siesta options
* ref@bsc--2007/siesta-bsc--master--2.1--patch-3
New geometry module
* ref@bsc--2007/siesta-bsc--master--2.1--patch-4
New stub module for sparse matrices
* ref@bsc--2007/siesta-bsc--master--2.1--patch-5
More options for running tests
* ref@bsc--2007/siesta-bsc--master--2.1--patch-6
Encapsulation of k-point handling
* ref@bsc--2007/siesta-bsc--master--2.1--patch-7
Initialize iza in struct_init
* ref@bsc--2007/siesta-bsc--master--2.1--patch-8
Siesta_todo slimming by M. Quero
* ref@bsc--2007/siesta-bsc--master--2.1--patch-9
Fix import of no_l in born_charge
* ref@bsc--2007/siesta-bsc--master--2.1--patch-10
Some conversions to Fortran90 by M. Quero
* ref@bsc--2007/siesta-bsc--master--2.1--patch-11
Fourth session at the BSC
* ref@bsc--2007/siesta-bsc--master--2.1--patch-12
Clarification of the scope of the stress variables
* ref@bsc--2007/siesta-bsc--master--2.1--patch-13
Creation of siesta_forces
* ref@bsc--2007/siesta-bsc--master--2.1--patch-14
Replacement of some allocatables by pointers and automatics
* ref@bsc--2007/siesta-bsc--master--2.1--patch-15
New m_energies and m_steps modules. Back to old k-point behavior
* ref@bsc--2007/siesta-bsc--master--2.1--patch-16
Merge of removal of integer and real variables from siesta_todo
* ref@bsc--2007/siesta-bsc--master--2.1--patch-17
Fixes for troublesome bugs in reference code
* ref@bsc--2007/siesta-bsc--master--2.1--patch-18
Final cleanup of siesta_todo
* ref@bsc--2007/siesta-bsc--master--2.1--patch-19
Add character(len=*) routines to alloc.F90
* ref@bsc--2007/siesta-bsc--master--2.1--patch-20
Pointers in fixed and setspatial. si2x1h test added to bsc-Makefile
* ref@bsc--2007/siesta-bsc--master--2.1--patch-21
Re-organization of pulay module
* ref@bsc--2007/siesta-bsc--master--2.1--patch-22
Reorganization of hsparse/xijorb calls with new neighbor module
* ref@bsc--2007/siesta-bsc--master--2.1--patch-23
VPATH-aware compilation for multiple executable versions
* ref@bsc--2007/siesta-bsc--master--2.1--patch-24
Allocatables to pointers I
* ref@bsc--2007/siesta-bsc--master--2.1--patch-25
Allocatables to pointers II
* ref@bsc--2007/siesta-bsc--master--2.1--patch-26
Allocatables to pointers III
* ref@bsc--2007/siesta-bsc--master--2.1--patch-27
Allocatables to pointers IV -- new neighbor code + nspecies fix
* ref@bsc--2007/siesta-bsc--master--2.1--patch-28
Explicit array-ness in calls in initatom and cellxc
* ref@bsc--2007/siesta-bsc--master--2.1--patch-29
Fix wrong allocations in cellxc.F
* ref@bsc--2007/siesta-bsc--master--2.1--patch-30
Explicit array extents in initatom.f
* ref@bsc--2007/siesta-bsc--master--2.1--patch-31
Avoid shrinking of density-matrix arrays for extrapol.
* ref@bsc--2007/siesta-bsc--master--2.1--patch-32
Fix typo in state_init.F
* ref@bsc--2007/siesta-bsc--master--2.1--patch-33
Execute SCF loop when nscf=1
* ref@bsc--2007/siesta-bsc--master--2.1--patch-34
Clarify bounds of SCF loop in siesta_forces.
* ref@bsc--2007/siesta-bsc--master--2.2--base-0
tag of ref@bsc--2007/siesta-bsc--master--2.1--patch-31
* ref@bsc--2007/siesta-bsc--master--2.3--base-0
tag of ref@bsc--2007/siesta-bsc--master--2.2--base-0
* ref@bsc--2007/siesta-bsc--master--2.3--patch-1
Prepare CHANGES file
* ref@bsc--2007/siesta-bsc--master--2.3--patch-2
Merge patch-log for patch-30 from 2.1
* ref@bsc--2007/siesta-bsc--master--2.3--patch-3
New treatment of fractional atoms in VCA. Bug fix in lmxo
* ref@bsc--2007/siesta-bsc--master--2.3--patch-4
Re-enabling of kgrid update in variable-cell calculations
* ref@bsc--2007/siesta-bsc--master--2.3--patch-5
Fix typo in state_init.F
* ref@bsc--2007/siesta-bsc--master--2.3--patch-6
Merge filtering package by Eduardo Anglada
* ref@bsc--2007/siesta-bsc--master--2.3--patch-7
Add graphite_c6_full test for more realistic vdW test
* ref@bsc--2007/siesta-bsc--master--2.3--patch-8
Add patchlog for k-point fix
* ref@bsc--2007/siesta-bsc--master--2.3--patch-9
Fix of zmatrix code to deal with degenerate case
* ref@bsc--2007/siesta-bsc--master--2.3--patch-10
Units conversion in Util/Optical/optical.f. Scripts. Cosmetics
* ref@bsc--2007/siesta-bsc--master--2.3--patch-11
Implementation of basis_enthalpy calculation. Zmatrix dependency
* ref@bsc--2007/siesta-bsc--master--2.3--patch-12
Option to use fractional rc's for multiple zeta
* ref@bsc--2007/siesta-bsc--master--2.3--patch-13
Sync vpath changes. Update siesta.tex
* ref@bsc--2007/siesta-bsc--master--2.3--patch-14
Fix tag in MPI send/receive in mulliken
* ref@bsc--2007/siesta-bsc--master--2.3--patch-15
Option to use fractional rc's for multiple zeta
* ref@bsc--2007/siesta-bsc--master--2.3--patch-16
Merge of XML tester by Eduardo Anglada. Portability fixes
* ref@bsc--2007/siesta-bsc--master--2.3--patch-17
Fix marenostrum-mpi.make. Syntax in compare_m.f90
* ref@bsc--2007/siesta-bsc--master--2.3--patch-18
MareNostrum fixes. zdrot to blas. obj_setup. compare_m syntax
* ref@bsc--2007/siesta-bsc--master--2.3--patch-19
Sync to bsc--2.1. Change banner in CHANGES to bsc--2.3
* ref@bsc--2007/siesta-bsc--master--2.3--patch-20
Execute SCF loop when nscf=1
* ref@bsc--2007/siesta-bsc--master--2.3--patch-21
Refinements of XML tester. New ioncat program
* ref@bsc--2007/siesta-bsc--master--2.3--patch-22
Undef var in kgridinit, iohs MPI write, pdosg array bound, Origin shift
* ref@bsc--2007/siesta-bsc--master--2.3--patch-23
Zmatrix optimization enhancements. Sign change in MM stress.
* ref@bsc--2007/siesta-bsc--master--2.3--patch-24
Clarify bounds of SCF loop in siesta_forces.
* ref@bsc--2007/siesta-bsc--master--2.3--patch-25
Avoid MP-grid permutations in trivial gamma case
* ref@bsc--2007/siesta-bsc--master--2.3--patch-26
Sync to reference--2.1
* siesta@uam.es--2006/siesta-bsc--reference--2.1--base-0
tag of siesta@uam.es--2006/siesta-devel--reference--2.1--patch-29
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-1
First stage of siesta.F splitting at BSC
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-2
Created struct_init for initial geometry setup. New test force_2
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-3
Consolidate geometry updates at the end of loop
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-4
Initialize vol2 correctly in m_check_supercell.f
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-5
Kgrid setup streamlined. Bands. Proximity check. Hsparse allocation
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-6
Work by Manuel Quero before the meeting on Feb 7th
* siesta@uam.es--2006/siesta-bsc--reference--2.1--patch-7
Moved the Born-effective-charge code
- files added:
- Docs/BSC-Issues
- Tests/force_2
- Tests/force_constants
- files removed:
- Src/Sys/intel-checks.make
- files renamed:
- Src/ksvinit.f => Src/ksvinit.F
- Src/pulayx.F => Src/m_pulay.F
- Src/zheevds.f => Src/zheevds.F90
- Tests/Makefile => Tests/std-Makefile
- Tests/test.mk => Tests/std-test.mk
- files modified:
- Docs/CHANGES
added
removed
Src/siesta.F
1
!
2
! This file is part of the SIESTA package.
3
!
4
! Copyright (c) Fundacion General Universidad Autonoma de Madrid:
5
! E.Artacho, J.Gale, A.Garcia, J.Junquera, P.Ordejon, D.Sanchez-Portal
6
! and J.M.Soler, 1996-2006.
7
!
8
! Use of this software constitutes agreement with the full conditions
9
! given in the SIESTA license, as signed by all legitimate users.
10
!
11
1
Program SIESTA
12
2
13
C ***
14
C SIESTA Density Functional LCAO program package.
15
C Copyright by Fundacion General Universidad Autonoma de Madrid:
16
C E.Artacho, J.Gale, A.Garcia, J.Junquera, P.Ordejon, D.Sanchez-Portal
17
C and J.M.Soler, 1996-2006
18
C ***
19
C Copy or disemination of all or part of this package is not permitted
20
C without prior and explicit authorization by the authors.
21
C Send comments/suggestions/bug-reports to siesta@uam.es
22
C ***
23
24
C Modules
25
use precision, only: dp
26
use parallel, only: Node, Nodes, IOnode, ParallelOverK
27
use m_cell, only: ucell
28
use atmfuncs, only: rcut, uion, izofis
29
use atomlist, only: xa, xalast, indxuo, na_u
30
use atomlist, only: isa, cisa, elem, indxua, iza, lastkb, lasto
31
use atomlist, only: qa, amass, qtot, no_u, rmaxo, rmaxkb
32
use atomlist, only: zvaltot
33
use atomlist, only: rmaxv, na_s, iphorb, iphkb, rco, rckb
34
use atomlist, only: datm, no_s, iaorb
35
use atomlist, only: superx, superc, initatomlists
36
!
37
use kpoint_grid, only: maxk, nkpnt, setup_kpoint_grid
38
use kpoint_grid, only: kpoint, kweight, writek
39
!
40
use band, only: initbands, bands
41
use fdf, only : fdf_block, fdf_convfac, fdf_boolean
42
use fdf, only : fdf_physical
43
use m_fdf_global, only: fdf_global_get
44
use sys, only: die, bye
45
use parse
46
use periodic_table, only : symbol
47
use xcmod, only: setXC
48
use molecularmechanics, only : inittwobody, twobody
49
use metaforce, only: lMetaForce, initmeta, meta
50
#ifdef MPI
51
use mpi_siesta, only: mpi_comm_world
52
use m_mpi_utils, only : globalize_sum, globalize_max
53
#endif
54
use m_mpi_utils, only : broadcast
55
use alloc, only: re_alloc, alloc_report
56
use phonon, only: phonon_num_disps, phonon_setup
57
use phonon, only: phonon_write_forces, phonon_restore_coords
58
use phonon, only: phonon_set_coords
59
60
use densematrix
61
62
use m_ordern, only: ordern
63
use m_hsparse, only: hsparse
64
65
use parallelsubs, only: getnodeorbs
66
use writewave, only: initwave, wwave
67
use iopipes, only: forcestopipe, coordsfrompipe
68
use m_iostruct, only: write_struct, read_struct
69
use m_nlefsm, only: nlefsm
70
use m_overfsm, only: overfsm
71
use m_check_supercell, only: check_sc_factors
72
73
use m_broyden_mixing
74
75
use siesta_cml
76
77
use files, only : slabel, label_length
78
79
use m_timestamp, only : timestamp
80
use m_wallclock, only : wallclock
81
use units
82
use zmatrix, only: lUseZmatrix, iofaZmat
83
use zmatrix, only: CartesianForce_to_ZmatForce
84
use zmatrix, only: write_Zmatrix,
85
$ write_canonical_ucell_and_Zmatrix
86
use m_broyden_optim, only:broyden_optimizer
87
use m_zm_broyden_optim, only: zm_broyden_optimizer
88
use m_cell_broyden_optim, only: cell_broyden_optimizer
89
use m_redata, only: redata
90
use m_ioxv, only: ioxv
91
92
use m_smearing, only: temp
93
use m_dynamics, only: nose, verlet2, npr, anneal, pr
94
use md_out, only: md_v_format
95
#ifdef CDF
96
use md_out, only: md_netcdf
97
#endif
98
use basis_enthalpy, only: write_basis_enthalpy
3
use m_siesta_init
4
use m_siesta_analysis
5
use m_siesta_move
6
use m_siesta_end
7
use m_siesta_forces
8
USE m_steps, only: inicoor, fincoor
99
9
100
10
implicit none
101
11
102
integer
103
. fincoor,
104
. i, ia, ia1, ia2, iadispl, ianneal,
105
. idyn, ifa, ifinal, ihmat, ihuge, iiscf, ik,
106
. ind, inicoor, io, ioptlwf, iord,
107
. iquench, is, isel, iscf,
108
. isolve, ispin, istp, istpsave, istart, istep, istr,
109
. iu, iunit, iv, ix, ixdispl,
110
. j, ja, jamin, jo, jx,
111
. level, maxbk, maxnh,
112
. maxna, maxsav, broyden_maxit, maxwk, mscell(3,3), mullipop,
113
. nkick, nnamax, nauxpul, nbcell, nbk, ncgmax, nh,
114
. nmove, nnia, ns, nsc(3), nscold(3), nscf, nspin,
115
. ntm(3), ntcon, no_l, nxij, nwk, pmax, nkpol, nhist,
116
. neigwanted, neigmin
117
118
#ifdef MPI
119
integer
120
. MPIerror, ntmp
121
#endif
122
123
integer, dimension(:), allocatable ::
124
. jna
125
126
integer, pointer, save ::
127
. listh(:), listhold(:), listhptr(:), listhptrold(:),
128
. numh(:), numhold(:)
129
130
real(dp)
131
. bcell(3,3), beta, bulkm,
132
. charnet, cfmax, cftem, const, cstress(3,3),
133
. mstress(3,3), kin_stress(3,3), ps(3,3),
134
. dipol(3), dDmax, dDtol, dEmax, dEtol, DEharr, DEna,
135
. dt, DUext, DUscf, Dxc, dx, dxmax, e1, e2, sigma,
136
. Ecorrec, ef, Eharrs, Eharrs1, Eions, Ekin, Ekinion, Emad, Ena,
137
. Elast, Enaatm, Enascf, Enl, Entrop, Entropy, Emeta,
138
. eta(2), etol, Etot, Etot_output, Exc, E0, Emm,
139
. factor, fmax, fmean, FreeE, fres, ftem, ftol, ftot(3),
140
. g2cut, g2max, getot, kn, kpr,
141
. mn, mpr, occtol, Pint, Pmol, Psol, rmax_bonds, Press
142
real(dp)
143
. qaux, qspin(4), qsol,
144
. rcoor, rcoorcp, rijmin, rmax, rmaxh, rmin, r2min,
145
. scell(3,3), stot, stress(3,3), strtol, svec(3),
146
. taurelax, tempinit, tempion, tiny, tp, ts, tt,
147
. tstress(3,3), Uatm, Uscf,
148
. vcell(3,3), virial, vn,
149
. volcel, volume, vpr, wmix, wmixkick,
150
. stressl(3,3), veclen
151
152
real(dp), dimension(:), allocatable ::
153
. efs, qs, r2ij
154
155
real(dp), pointer, save ::
156
. H0(:), S(:), wgthpol(:)
157
158
real(dp), dimension(:,:), allocatable ::
159
. auxpul, cfa, fa, fal,
160
. polR, polxyz, va, xij, xijo, wfk
161
162
real(dp), pointer, save ::
163
. Dold(:,:), Dscf(:,:), Dscfsave(:,:), Eold(:,:), Escf(:,:),
164
. bk(:,:), H(:,:), kpol(:,:)
165
166
#ifdef MPI
167
real(dp)
168
. buffer1, stresstmp(3,3), qtmp(4)
169
170
real(dp), dimension(:,:), allocatable ::
171
. fatmp
172
#endif
173
174
real(dp), dimension(:,:,:), allocatable :: ebk
175
real(dp), dimension(:,:,:), pointer :: eo, qo
176
177
logical
178
. auxchanged, chebef, dminit, default, dumpcharge,
179
. final, first, fixauxcell, fixspin, found, foundxv, foundzm,
180
. gamma, naiveauxcell,
181
. initdmaux, inspn, itest, last, lastst, mix, mmix, negl, noeta,
182
. outlng, overflow, overflowed, pulfile, relaxd,
183
. savehs, savevh, savevt, savdrh, savrho,
184
. savepsch, savetoch,
185
. usesavecg, usesavelwf, usesavedm, usesavedmloc, usesavexv,
186
. usesavezm, writeig,
187
. writbk, writmd, writpx, writb, writec, writef,
188
. writic, varcel, genlogic, do_pdos, writedm, atmonly,
189
. harrisfun, muldeb, eggbox_block, use_struct_file,
190
$ require_energy_convergence, broyden_optim, struct_only,
191
$ bornz, cell_can_change, change_kgrid_in_md,
192
$ RelaxCellOnly, RemoveIntraMolecularPressure
193
194
character
195
. line*150, sname*150, shape*10, message*79
196
197
character(len=label_length+5) :: fildrh
198
character(len=label_length+5) :: filepsch
199
character(len=label_length+5) :: filetoch
200
character(len=label_length+5) :: filevh
201
character(len=label_length+5) :: filevt
202
character(len=label_length+5) :: filrho
203
character(len=label_length+5) :: paste
204
205
external :: automatic_cell,
206
. bonds, cgvc, cgvc_zmatrix, fixed,
207
. dhscf, diagon, dnaefs, extrapol, initatom,
208
. iodm, iozm,
209
. kinefsm, mulliken, naefs, neighb,
210
. paste, pulayx,
211
. reinit, shaper, spnvec,
212
. timer, volcel, xijorb, memory,
213
. ioeig, iofa, iokp, iomd, prversion, eggbox
214
215
type(parsed_line), pointer :: p
216
217
data
218
. e1, e2 / 1._dp, -1._dp /
219
. eggbox_block /.true./
220
. final /.false./
221
. relaxd /.false./
222
. tiny /1.e-15_dp/
223
. ihuge /1073741823/
224
! maxnh is (initially) the minimum size of all the sparse arrays.
225
! It must be 1 (not 0) since we frequently use the F77 idiom of
226
! passing the first element of the array instead of all of it.
227
. maxnh / 1 /
228
. maxna / 200 /
229
. tempion / 0.0_dp /
230
. no_l, nnamax /2*1/
231
. nauxpul, nbk, ns, nspin, nxij /5*1/
232
. nscold / 3*0 /
233
. nsc, mscell / 1,1,1, 1,0,0, 0,1,0, 0,0,1 /
234
c---------------------------------------------------------------------
235
236
C Initialise MPI and set processor number
237
#ifdef MPI
238
call MPI_Init( MPIerror )
239
call MPI_Comm_Rank( MPI_Comm_World, Node, MPIerror )
240
call MPI_Comm_Size( MPI_Comm_World, Nodes, MPIerror )
241
#endif
242
243
IOnode = (Node .eq. 0)
244
245
C Print version information ...........................................
246
if (IOnode) then
247
call prversion
248
#ifdef MPI
249
if (Nodes.gt.1) then
250
write(6,'(/,a,i4,a)')
251
. '* Running on ', Nodes, ' nodes in parallel'
252
else
253
write(6,'(/,a,i4,a)')
254
. '* Running in serial mode with MPI'
255
endif
256
#else
257
write(6,'(/,a,i4,a)')
258
. '* Running in serial mode'
259
#endif
260
call timestamp('Start of run')
261
call wallclock('Start of run')
262
endif
263
C ..................
264
265
C Start time counter ..................................................
266
call timer( 'siesta', 0 )
267
call timer( 'siesta', 1 )
268
call timer( 'Setup', 1 )
269
270
C Nullify arrays
271
nullify(Haux,Saux,psi)
272
nullify(bk,H,H0,S,kpol,wgthpol)
273
nullify(Dold,Dscf,Dscfsave,Eold,Escf)
274
nullify(listh,listhold,listhptr,listhptrold)
275
nullify(numh,numhold)
276
277
C Initialize some variables
278
DUext = 0.0_dp
279
Eharrs = 0.0_dp
280
Eharrs1 = 0.0_dp
281
Eions = 0.0_dp
282
Ekinion = 0.0_dp
283
Elast = 0.0_dp
284
Emad = 0.0_dp
285
Emeta = 0.0_dp
286
Emm = 0.0_dp
287
Entrop = 0.0_dp
288
Entropy = 0.0_dp
289
FreeE = 0.0_dp
290
291
C Initialise read .....................................................
292
call reinit(sname)
293
294
call siesta_cml_init() ! Initialize CML (relies on reinit)
295
296
C Set allocation report level .........................................
297
call fdf_global_get(level, 'alloc_report_level', 0 )
298
call alloc_report( level=level, file=trim(slabel)//'.alloc' )
299
C ..................
300
301
C Initialise exchange-correlation functional information
302
call setXC()
303
304
C Initialise force field component
305
call inittwobody
306
307
C Initialize pseudopotentials and atomic orbitals
308
if (IOnode) call initatom()
309
call broadcast_basis()
310
311
call fdf_global_get(atmonly,'Atom-Setup-Only',.false.)
312
if (atmonly) call bye("End of atom setup")
313
314
if (Node.eq.0) then
315
write(6,'(/,a,20("*"),a,28("*"))')
316
. 'siesta: ', ' Simulation parameters '
317
write(6,'(a)') 'siesta:'
318
write(6,'(a)') 'siesta: The following are some of the '//
319
. 'parameters of the simulation.'
320
write(6,'(a)') 'siesta: A complete list of the parameters '//
321
. 'used, including default values,'
322
write(6,'(a,a)')'siesta: can be found in file out.fdf'
323
write(6,'(a)') 'siesta:'
324
endif
325
326
C Read simulation sizes ...............................................
327
328
!! Read number of atoms and coordinates, and unit cell
329
330
call fdf_global_get(use_struct_file,"MD.UseStructFile",.false.)
331
if (use_struct_file) then
332
call read_struct( na_u, ucell) ! Sets na_u, xa, and isa
333
else
334
call coor(na_u,ucell) ! Sets na_u, xa, and isa
335
endif
336
337
338
if (cml_p) then
339
C We need the names of the elements on node 0
340
allocate(elem(na_u))
341
call memory('A','S',2*na_u,'siesta')
342
do i = 1, na_u
343
elem(i) = symbol(izofis(isa(i)))
344
enddo
345
endif
346
347
C Allocate arrays based on read sizes ................
348
!
349
! Those living in module atomlists
350
!
351
nullify(indxua,iza,lastkb,lasto,qa,amass,xalast)
352
call re_alloc(indxua,1,na_u,name='indxua',routine='siesta')
353
call re_alloc(iza,1,na_u,name='iza',routine='siesta')
354
call re_alloc(lastkb,0,na_u,name='lastkb',routine='siesta')
355
call re_alloc(lasto,0,na_u,name='lasto',routine='siesta')
356
call re_alloc(qa,1,na_u,name='qa',routine='siesta')
357
call re_alloc(xalast,1,3,1,na_u,name='xalast',routine='siesta')
358
call re_alloc(amass,1,na_u,name='amass',routine='siesta')
359
360
!
361
! Others
362
!
363
allocate(cfa(3,na_u))
364
call memory('A','D',3*na_u,'siesta')
365
allocate(fa(3,na_u))
366
call memory('A','D',3*na_u,'siesta')
367
allocate(fal(3,na_u))
368
call memory('A','D',3*na_u,'siesta')
369
allocate(va(3,na_u))
370
call memory('A','D',3*na_u,'siesta')
371
372
C Initialise those arrays that must be pre-initialised
373
fal(1:3,1:na_u) = 0.0_dp
374
stress(1:3,1:3) = 0.0_dp
375
376
if (IOnode) call spin_init(nspin)
377
call broadcast(nspin)
378
379
allocate(efs(nspin))
380
call memory('A','D',nspin,'siesta')
381
allocate(qs(nspin))
382
call memory('A','D',nspin,'siesta')
383
allocate(polR(3,nspin))
384
call memory('A','D',3*nspin,'siesta')
385
allocate(polxyz(3,nspin))
386
call memory('A','D',3*nspin,'siesta')
387
388
C Read simulation data ................................................
389
call redata(na_u, ns, nspin, outlng, g2cut, charnet, negl, nscf,
390
. dDtol, dEtol, mix, wmix, isolve, temp, fixspin, ts,
391
. ncgmax, ftol, strtol, eta, etol, rcoor,
392
. ioptlwf, chebef, noeta, rcoorcp, beta, pmax,
393
. idyn, istart, ifinal, nmove, ianneal, iquench,
394
. dt, ia1, ia2, dx, dxmax, tt, tp, mn, mpr,
395
. bulkm, taurelax,
396
. usesavelwf, usesavedm, usesavecg,
397
. mullipop, inspn, maxsav, nkick, wmixkick,
398
. pulfile, tempinit, dumpcharge, varcel, harrisfun,
399
. occtol,broyden_maxit,require_energy_convergence,
400
$ broyden_optim)
401
402
C Find some switches ..................................................
403
call fdf_global_get(writek,'WriteKpoints' , outlng )
404
call fdf_global_get(writef,'WriteForces' , outlng )
405
call fdf_global_get(writedm,'WriteDM' , .true.)
406
call fdf_global_get(writb, 'WriteBands' , outlng )
407
call fdf_global_get(writbk, 'WriteKbands' , outlng )
408
call fdf_global_get(writeig,'WriteEigenvalues', outlng )
409
call fdf_global_get(writec, 'WriteCoorStep' , outlng )
410
call fdf_global_get(writic, 'WriteCoorInitial', .true. )
411
call fdf_global_get(writmd, 'WriteMDhistory' , .false.)
412
call fdf_global_get(writpx, 'WriteMDXmol' , .not. writec)
413
call fdf_global_get(default, 'UseSaveData' , .false.)
414
call fdf_global_get(usesavexv, 'MD.UseSaveXV' , default)
415
call fdf_global_get(usesavezm, 'MD.UseSaveZM' , default)
416
call fdf_global_get(savehs, 'SaveHS' , .false.)
417
call fdf_global_get(fixauxcell, 'FixAuxiliaryCell', .false.)
418
call fdf_global_get(naiveauxcell, 'NaiveAuxiliaryCell', .false.)
419
call fdf_global_get(initdmaux, 'ReInitialiseDM' , .false.)
420
call fdf_global_get(muldeb, 'MullikenInSCF' , .false.)
421
call fdf_global_get(rijmin, 'WarningMinimumAtomicDistance',
422
. 1.0_dp, 'Bohr' )
423
call fdf_global_get(bornz, 'BornCharge' , .false.)
424
if (idyn.ne.6) bornz = .false.
425
426
call fdf_global_get(change_kgrid_in_md,
427
$ "ChangeKgridInMD", .false.)
428
call fdf_global_get(ParallelOverK, 'Diag.ParallelOverK', .false.)
429
call fdf_global_get(RelaxCellOnly, 'MD.RelaxCellOnly', .false.)
430
call fdf_global_get(RemoveIntraMolecularPressure,
431
$ 'MD.RemoveIntraMolecularPressure', .false.)
432
433
C Read Z-matrix coordinates and forces from file
434
if (lUseZmatrix) then
435
foundzm = .false.
436
if (usesavezm) then
437
call iozm('read',ucell,vcell,xa,foundzm)
438
if (IOnode) then
439
if (.not.foundzm) then
440
write(6,'(/,a)') 'siesta: WARNING: ZM file not found'
441
else
442
write(6,'(/,a)')
443
. "! Info in XV file prevails over previous structure input"
444
endif
445
endif
446
endif
447
endif
448
449
C Read cell shape and atomic positions from a former run ..............
450
foundxv = .false.
451
if (usesavexv) then
452
call ioxv('read', ucell, vcell, na_u, isa, iza, xa, va, foundxv)
453
if (IOnode) then
454
if (.not.foundxv) then
455
write(6,'(/,a)') 'siesta: WARNING: XV file not found'
456
else
457
write(6,"(a)")
458
$ "! Info in XV file prevails over previous structure input"
459
endif
460
endif
461
endif
462
C ..................
463
464
C Read cell shape and atomic positions from driver program through pipe
465
if (idyn.eq.8) then
466
call coordsFromPipe( na_u, xa, ucell )
467
end if
468
C .....................
469
470
C Dump initial coordinates to output ..................................
471
if ( writic.and.(IOnode) ) then
472
write(6,'(/a)') 'siesta: Atomic coordinates (Bohr) and species'
473
write(6,"('siesta: ',2x,3f10.5,i3,3x,i6)")
474
. ( (xa(ix,ia), ix=1,3), isa(ia), ia, ia=1, na_u)
475
endif
476
C ..................
477
478
479
C Initialize atom lists
480
call initatomlists() ! Sets iza
481
qtot = qtot - charnet
482
if (IOnode) then
483
write(6,fmt="(a,f12.6)") "Total number of electrons: ", qtot
484
write(6,fmt="(a,f12.6)") "Total ionic charge: ", zvaltot
485
endif
486
487
488
C Calculate spin populations for fixed spin case...
489
if (fixspin) then
490
if (nspin .ne. 2)
491
$ call die('siesta: ERROR: ' //
492
$ 'You can only fix the spin of the system' //
493
$ ' for collinear spin polarized calculations.')
494
do i = 1,2
495
qs(i) = (qtot + (3-2*i)*ts) / 2.0_dp
496
enddo
497
else
498
qs(1:nspin) = 0.0_dp
499
if (nspin .le. 2) then
500
do ispin = 1,nspin
501
qs(ispin) = qtot/nspin
502
enddo
503
endif
504
endif
505
C ..................
506
507
C Find maximum interaction range ......................................
508
if (negl) then
509
rmaxh = 2.0_dp*rmaxo
510
else
511
rmaxh = 2.0_dp*rmaxo + 2.0_dp*rmaxkb
512
endif
513
C ......................
514
515
C Automatic cell generation ...........................................
516
if (volcel(ucell) .lt. 1.0d-8) then
517
call automatic_cell(ucell,scell,na_u,xa,isa,charnet)
518
endif
519
520
C Find system shape ...................................................
521
call shaper( ucell, na_u, isa, xa, shape, nbcell, bcell )
522
if (IOnode) then
523
write(6,'(/,2a)') 'siesta: System type = ', shape
524
endif
525
526
C Find interatomic distances (output in file BONDS)
527
if (IOnode) then
528
rmax_bonds = fdf_physical("MaxBondDistance", 6.0_dp, "Bohr")
529
call bonds( ucell, na_u, isa, xa,
530
$ rmax_bonds, trim(slabel)// ".BONDS" )
531
endif
532
533
C Output of initial system details:
534
if (cml_p) then
535
call cmlStartModule(xf=mainXML, title='Initial System')
536
call cmlAddMolecule(xf=mainXML, natoms=na_u,
537
. coords=xa/Ang, elements=elem, refs=cisa,
538
. style='x3', fmt='(f12.6)')
539
call cmlAddLattice(xf=mainXML, cell=ucell,
540
. units='siestaUnits:angstrom', dictref='siesta:ucell')
541
call cmlAddProperty(xf=mainXML, property=trim(shape),
542
. dictref='siesta:shape')
543
call cmlEndModule(xf=mainXML)
544
endif
545
546
! early exit if only checking the structure
547
548
call fdf_global_get(struct_only,'Output-Structure-Only',.false.)
549
if (IONode) then
550
call write_struct( ucell, na_u, isa, iza, xa )
551
if (lUseZmatrix) then
552
call write_canonical_ucell_and_Zmatrix()
553
endif
554
endif
555
if (struct_only) then
556
call bye("End of structure processing")
557
endif
558
559
!!! --------------------- End of Structure Generation section
560
561
C Madelung correction for charged systems .............................
562
if (charnet .ne. 0.0_dp) then
563
call madelung(ucell, shape, charnet, Emad)
564
endif
565
566
C Parallel initialisation
567
call initparallel(no_u,na_u,lasto,xa,ucell,rmaxh,rcoor,isolve)
568
if (IOnode) call show_distribution()
569
570
C Find number of locally stored orbitals and allocated related arrays
571
call GetNodeOrbs(no_u,Node,Nodes,no_l)
572
573
C Initialise arrays
574
call re_alloc(listhptr,1,no_l,name='listhptr',routine='siesta',
575
. copy=.false.)
576
call re_alloc(listhptrold,1,no_l,name='listhptrold',
577
. routine='siesta',copy=.false.)
578
call re_alloc(numh,1,no_l,name='numh',routine='siesta',
579
. copy=.false.)
580
call re_alloc(numhold,1,no_l,name='numhold',routine='siesta',
581
. copy=.false.)
582
listhptr(1:no_l) = 0
583
listhptrold(:) = 0
584
numh(:) = 0
585
numhold(:) = 0
586
587
C Get number of eigenstates that need to be calculated
588
589
call fdf_global_get(neigwanted,'NumberOfEigenStates',no_u)
590
591
592
C Check number of eigenstates - cannot be larger than number of
593
C basis functions or smaller than number of occupied states + 1
594
C so that the Fermi level can be estimated
595
do is = 1,nspin
596
neigmin = nint(qs(is)/real(3 - min(nspin,2), kind=dp)) + 1
597
neigwanted = max(neigwanted,neigmin)
598
enddo
599
neigwanted = min(neigwanted,no_u)
600
601
! Find k-grid for Brillouin zone integration
602
! Sets kscell, nkpnt, maxk, kpoint, kweight
603
call setup_Kpoint_grid( ucell )
604
605
nullify(eo,qo)
606
call re_alloc(eo,1,no_u,1,nspin,1,maxk,name="eo",
607
$ routine="state_init")
608
call re_alloc(qo,1,no_u,1,nspin,1,maxk,name="qo",
609
$ routine="state_init")
610
611
C ......................
612
613
C Find number of band k-points ........................................
614
nbk = 0
615
maxbk = 1
616
call re_alloc(bk,1,3,1,maxbk,name='bk',routine='siesta',
617
. copy=.false.)
618
C
619
call initbands( maxbk, nbk, bk )
620
C
621
if (nbk .gt. maxbk) then
622
C If there wasn't enough space to store bands on first call correct
623
C the dimensions and repeat the initialisation
624
maxbk = max(nbk,1)
625
call re_alloc(bk,1,3,1,maxbk,name='bk',routine='siesta',
626
. copy=.false.)
627
nbk = 0
628
call initbands( maxbk, nbk, bk )
629
endif
630
allocate(ebk(no_u,nspin,maxbk))
631
call memory('A','D',no_u*nspin*maxbk,'siesta')
632
633
C ......................
634
635
C Find number of k-points for wavefunction printout ....................
636
nwk = 0
637
maxwk = 1
638
allocate(wfk(3,maxwk))
639
call memory('A','D',3*maxwk,'siesta')
640
641
call initwave( maxwk, no_u, nwk, wfk, overflow )
642
643
if (overflow) then
644
if (nwk .gt. maxwk) then
645
C If there wasn't enough space to store bands on first call correct
646
C the dimensions and repeat the initialisation
647
maxwk = max(nwk,1)
648
call memory('D','D',size(wfk),'siesta')
649
deallocate(wfk)
650
allocate(wfk(3,maxwk))
651
call memory('A','D',3*maxwk,'siesta')
652
endif
653
nwk = 0
654
call initwave( maxwk, no_u, nwk, wfk, overflow )
655
if (overflow)
656
$ call die('siesta: ERROR: Unsuccessful initialization of' //
657
$ ' list of wavefunctions to print')
658
endif
659
C ......................
660
661
C Find the grid for the calculation of the polarization..............
662
nkpol = 1
663
call re_alloc(kpol,1,3,1,nkpol,name='kpol',routine='siesta',
664
. copy=.false.)
665
call re_alloc(wgthpol,1,nkpol,name='wgthpol',routine='siesta',
666
. copy=.false.)
667
668
call KSV_init(ucell, 0, nkpol, kpol, wgthpol)
669
670
call re_alloc(kpol,1,3,1,nkpol,name='kpol',routine='siesta',
671
. shrink=.false.,copy=.false.)
672
call re_alloc(wgthpol,1,nkpol,name='wgthpol',routine='siesta',
673
. shrink=.false.,copy=.false.)
674
675
C Find if only gamma point is used ....................................
676
if (nkpnt.eq.1 .and. abs(kpoint(1,1)).lt.tiny .and.
677
. abs(kpoint(2,1)).lt.tiny .and.
678
. abs(kpoint(3,1)).lt.tiny) then
679
gamma = .true.
680
else
681
gamma = .false.
682
endif
683
if (nbk .gt. 0) gamma = .false.
684
if (nwk .gt. 1) gamma = .false.
685
if (nwk .eq. 1) then
686
if (abs(wfk(1,1)).gt.tiny .and.
687
. abs(wfk(2,1)).gt.tiny .and.
688
. abs(wfk(3,1)).gt.tiny) then
689
gamma = .false.
690
endif
691
endif
692
if (nkpol.gt.0) gamma = .false.
693
C ....................
694
695
C Find required supercell
696
C 2*rmaxh is used to guarantee that two given orbitals in the
697
C supercell can only overlap once
698
699
if (gamma) then
700
nsc(1:3) = 1
701
else
702
do i=1,3
703
veclen = sqrt(ucell(1,i)**2+ucell(2,i)**2+ucell(3,i)**2)
704
nsc(i) = ceiling( 2 * rmaxh / veclen )
705
enddo
706
if (.not. naiveauxcell)
707
$ call check_sc_factors(ucell,nsc,2*rmaxh)
708
endif
709
710
mscell = 0.0_dp
711
do i = 1, 3
712
mscell(i,i) = nsc(i)
713
nscold(i) = nsc(i)
714
enddo
715
716
717
C Find auxiliary supercell (required only for k sampling) ............
718
call superc( ucell, scell, nsc)
719
720
C Initialise metadynamic forces if required
721
call initmeta
722
723
C Initialize atomic velocities to zero ................................
724
if (.not. foundxv) then
725
va(1:3,1:na_u) = 0.0_dp
726
vcell(1:3,1:3) = 0.0_dp
727
endif
728
C ..................
729
730
C Begin of coordinate relaxation iteration ============================
731
C Notice that this loop is not indented
732
if (idyn .eq. 0) then
733
inicoor = 0
734
fincoor = nmove
735
else if (idyn .ge. 1 .and. idyn .le. 5) then
736
inicoor = istart
737
fincoor = ifinal
738
else if (idyn .eq. 6) then
739
inicoor = 0
740
fincoor = (ia2-ia1+1)*3*2
741
else if (idyn .eq. 7) then
742
call phonon_setup
743
inicoor = 1
744
fincoor = phonon_num_disps
745
else if (idyn .eq. 8) then
746
inicoor = 0
747
fincoor = ihuge
748
else
749
call die('siesta: wrong idyn')
750
endif
751
752
C Build initial velocities according to Maxwell-Bolzmann distribution....
753
if (idyn .ne. 0 .and. idyn .ne. 6 .and. (.not. foundxv))
754
. call vmb(na_u,tempinit,amass,xa,isa,va)
755
C ..................
756
757
istp = 0
758
call timer( 'Setup', 2 )
759
760
C Output memory use before main loop
761
call printmemory( 6, 0 )
762
763
C Initialization now complete. Flush stdout.
764
if (ionode) call pxfflush(6)
765
766
C Start loop over coordinate changes
767
768
do istep = inicoor,fincoor
769
call timer( 'IterMD', 1 )
770
istp = istp + 1
771
if (IOnode) then
772
write(6,'(/2a)') 'siesta: ',
773
. '=============================='
774
select case (idyn)
775
case (0)
776
write(6,'(28(" "),a,i6)') 'Begin CG move = ',istep
777
if (cml_p) call cmlStartStep(mainXML, type='CG', index=istp)
778
case (1:5)
779
write(6,'(28(" "),a,i6)') 'Begin MD step = ',istep
780
if (cml_p) call cmlStartStep(mainXML, type='MD', index=istep)
781
case (6)
782
write(6,'(28(" "),a,i6)') 'Begin FC step = ',istep
783
if (cml_p) call cmlStartStep(mainXML, type='FC', index=istep)
784
if (istep .eq. 0) then
785
write(6,'(28(" "),a)') 'Undisplaced coordinates'
786
else
787
iadispl = (istep-mod(istep-1,6))/6+ia1
788
write(6,'(28(" "),a,i6)') 'displace atom ',
789
. iadispl
790
ix = mod(istep-1,6)+1
791
ixdispl = (ix - mod(ix-1,2) +1)/2
792
write(6,'(28(" "),a,i6)') 'in direction ',
793
. ixdispl
794
dx=-dx
795
write(6,'(28(" "),a,f8.4,a)') 'by ',
796
. dx, ' Bohr'
797
C Displace atom by dx...
798
xa(ixdispl,iadispl) = xa(ixdispl,iadispl) + dx
799
endif
800
case (7)
801
call phonon_set_coords(istep,xa,ucell)
802
case (8)
803
write(6,'(28(" "),a,i6)') 'Begin Server step = ',istep
804
if (cml_p) call cmlStartStep(mainXML, type='FS', index=istep)
805
806
end select
807
write(6,'(2a)') ' ',
808
. '=============================='
809
else ! not IOnode
810
select case (idyn)
811
C We don't need to do anything for 0<idyn<6
812
case(6)
813
if (istep .ne. 0) then
814
iadispl = (istep-mod(istep-1,6))/6 + ia1
815
ix = mod(istep-1,6) + 1
816
ixdispl = (ix - mod(ix-1,2) +1)/2
817
dx = - dx
818
C Displace atom by dx...
819
xa(ixdispl,iadispl) = xa(ixdispl,iadispl) + dx
820
endif
821
case(7)
822
call phonon_set_coords(istep,xa,ucell)
823
end select
824
endif
825
826
C Get coordinates from driver program through pipe
827
if (idyn.eq.8 .and. istep.ne.inicoor) then
828
call coordsFromPipe( na_u, xa, ucell )
829
if (volcel(ucell) < 1.0e-8_dp) then
830
call automatic_cell(ucell,scell,na_u,xa,isa,charnet)
831
endif
832
end if
833
834
if (IOnode) then
835
C Print Z-matrix coordinates
836
if (lUseZmatrix) then
837
call write_Zmatrix
838
endif
839
C Print atomic coordinates ............................................
840
call outcoor( ucell, xa, na_u, ' ', writec )
841
call siesta_write_positions()
842
endif
843
C ...................
844
845
C Actualize things if unit cell might have changed
846
auxchanged = .false.
847
cell_can_change = ( varcel .or.
848
$ (idyn .eq. 8) ! Force/stress evaluation
849
$ )
850
if (change_kgrid_in_md) then
851
cell_can_change = cell_can_change .or.
852
$ (idyn .eq. 3) ! Parrinello-Rahman
853
$ .or. (idyn .eq. 4) ! Nose-Parrinello-Rahman
854
$ .or. (idyn .eq. 5) ! Anneal
855
endif
856
857
if ( cell_can_change
858
$ .and. (istep.ne.inicoor) .and. (.not.gamma) ) then
859
860
call setup_Kpoint_grid( ucell )
861
862
call re_alloc(eo,1,no_u,1,nspin,1,maxk,name="eo",
863
$ routine="state_init")
864
call re_alloc(qo,1,no_u,1,nspin,1,maxk,name="qo",
865
$ routine="state_init")
866
867
C Find required supercell
868
869
auxchanged = .false.
870
871
if (gamma) then
872
nsc(1:3) = 1
873
else if (fixauxcell) then
874
nsc(i) = nscold(i)
875
else
876
do i=1,3
877
veclen = sqrt(ucell(1,i)**2+ucell(2,i)**2+ucell(3,i)**2)
878
nsc(i) = ceiling( 2 * rmaxh / veclen )
879
enddo
880
if (.not. naiveauxcell)
881
$ call check_sc_factors(ucell,nsc,2*rmaxh)
882
endif
883
884
mscell = 0.0_dp
885
do i = 1, 3
886
mscell(i,i) = nsc(i)
887
if (nsc(i).ne.nscold(i)) auxchanged = .true.
888
nscold(i) = nsc(i)
889
enddo
890
891
C Madelung correction for charged systems .............................
892
if (charnet .ne. 0.0_dp) then
893
call madelung(ucell, shape, charnet, Emad)
894
endif
895
896
endif
897
C End variable cell actualization
898
899
C Auxiliary supercell
900
call superc(ucell, scell, nsc)
901
902
C Print unit cell and find its volume
903
if (IOnode) call outcell(ucell)
904
volume = volcel( ucell )
905
C ...................
906
907
C Initialize neighb subroutine ........................................
908
144 ia = 0
909
isel = 0
910
rmax = max( 2._dp*rmaxv, 2._dp*rmaxo, rmaxo+rmaxkb )
911
nnia = maxna
912
if (allocated(jna)) then
913
call memory('D','I',size(jna),'siesta')
914
deallocate(jna)
915
endif
916
if (allocated(r2ij)) then
917
call memory('D','D',size(r2ij),'siesta')
918
deallocate(r2ij)
919
endif
920
if (allocated(xij)) then
921
call memory('D','D',size(xij),'siesta')
922
deallocate(xij)
923
endif
924
allocate(jna(maxna))
925
call memory('A','I',maxna,'siesta')
926
allocate(r2ij(maxna))
927
call memory('A','D',maxna,'siesta')
928
allocate(xij(3,maxna))
929
call memory('A','D',3*maxna,'siesta')
930
call neighb( scell, rmax, na_s, xa, ia, isel,
931
. nnia, jna, xij, r2ij )
932
nnamax = 0
933
do ia = 1,na_s
934
nnia = 0
935
call neighb( scell, rmax, na_s, xa, ia, isel,
936
. nnia, jna, xij, r2ij )
937
nnamax = max( nnamax, nnia )
938
enddo
939
if (nnamax .gt. maxna) then
940
C Increase maxna with safety margin when atoms move
941
maxna = nnamax + 0.10 * nnamax + 10
942
overflow = .true.
943
else
944
overflow = .false.
945
endif
946
if (overflow) goto 144
947
C ..................
948
949
C Check if any two atoms are unreasonably close .......................
950
do ia = 1,na_s
951
r2min = huge(1._dp)
952
jamin = 0
953
nnia = maxna
954
call neighb( scell, rmax, na_s, xa, ia, isel,
955
. nnia, jna, xij, r2ij )
956
do j = 1,nnia
957
ja = jna(j)
958
if ( r2ij(j).lt.r2min .and. ja.ge.ia ) then
959
C Check that it is not the same atom
960
if ( ja.ne.ia .or. r2ij(j).gt.1.d-12 ) then
961
r2min = r2ij(j)
962
jamin = ja
963
endif
964
endif
965
enddo
966
rmin = sqrt( r2min )
967
if (IOnode) then
968
if ( rmin .lt. rijmin ) write(6,'(a,2i6,a,f12.6,a)')
969
. 'siesta: WARNING: Atoms', ia, jamin, ' too close: rij =',
970
. rmin/Ang, ' Ang'
971
endif
972
enddo
973
C ..................
974
975
C List of nonzero Hamiltonian matrix elements .........................
976
overflow=.true.
977
overflowed=.false.
978
do while (overflow)
979
nh = maxnh
980
call re_alloc(listh,1,maxnh,name='listh',routine='siesta',
981
. copy=.false.)
982
call hsparse( negl, scell, nsc, na_s, isa, xa, lasto, lastkb,
983
. iphorb, iphKB, nh, numh, listhptr, listh )
984
if (nh .gt. maxnh) then
985
! Increase maxnh with safety margin for when atoms move
986
maxnh = 1.05 * nh + 40
987
overflowed=.true.
988
else
989
overflow=.false.
990
endif
991
enddo
992
! In first step, allocate anyway (to catch corner case
993
! where one node has nh=0, and doesn't overflow)
994
if (istp==1) then
995
call re_alloc(Dscf,1,maxnh,1,nspin,name='Dscf',
996
. routine='siesta',copy=.false.)
997
call re_alloc(Dscfsave,1,maxnh,1,nspin,name='Dscfsave',
998
. routine='siesta',copy=.false.)
999
call re_alloc(listhold,1,maxnh,name='listhold',
1000
. routine='siesta',copy=.false.)
1001
! Initialise Dscfsave to avoid problems in extrapol
1002
Dscfsave(1:maxnh,1:nspin) = 0.0_dp
1003
listhold(1:maxnh) = 0
1004
elseif (overflowed) then
1005
! We need to preserve the contents of these arrays
1006
! when reallocating.
1007
call re_alloc(Dscf,1,maxnh,1,nspin,name='Dscf',
1008
. routine='siesta',copy=.true.)
1009
call re_alloc(Dscfsave,1,maxnh,1,nspin,name='Dscfsave',
1010
. routine='siesta',copy=.true.)
1011
call re_alloc(listhold,1,maxnh,name='listhold',
1012
. routine='siesta',copy=.true.)
1013
endif
1014
if (istp==1.or.overflowed) then
1015
call re_alloc(Dold,1,maxnh,1,nspin,name='Dold',
1016
. routine='siesta',copy=.false.)
1017
call re_alloc(Eold,1,maxnh,1,nspin,name='Eold',
1018
. routine='siesta',copy=.false.)
1019
call re_alloc(Escf,1,maxnh,1,nspin,name='Escf',
1020
. routine='siesta',copy=.false.)
1021
endif
1022
1023
C Allocate/reallocate storage associated with Hamiltonian/Overlap matrix
1024
call re_alloc(H,1,maxnh,1,nspin,name='H',
1025
. routine='siesta',shrink=.false.,copy=.false.)
1026
call re_alloc(H0,1,maxnh,name='H0',routine='siesta',
1027
. shrink=.false.,copy=.false.)
1028
call re_alloc(S,1,maxnh,name='S',routine='siesta',
1029
. shrink=.false.,copy=.false.)
1030
1031
C ..................
1032
1033
C Some printout for debugging ........................................
1034
* if (IOnode) then
1035
* write(6,'(/,a)') 'siesta: connected orbitals'
1036
* do io = 1,no_u
1037
* call GlobalToLocalOrb(io,Node,Nodes,iio)
1038
* if (iio.gt.0) then
1039
* write(6,'(i6,4x,15i4)')
1040
* . io, (listh(listhptr(iio)+j),j=1,numh(iio))
1041
* endif
1042
*#ifdef MPI
1043
* call MPI_Barrier(MPI_Comm_World,MPIerror)
1044
*#endif
1045
* enddo
1046
* write(6,*) ' '
1047
* endif
1048
C ..................
1049
1050
C Find vectors between orbital centers ................................
1051
if (allocated(xijo)) then
1052
call memory('D','D',size(xijo),'siesta')
1053
deallocate(xijo)
1054
endif
1055
if (.not.gamma) then
1056
nxij = maxnh
1057
allocate(xijo(3,nxij))
1058
call memory('A','D',3*nxij,'siesta')
1059
call xijorb( negl, scell, na_u, na_s, xa,
1060
. lasto, lastkb, rco, rckb,
1061
. maxnh, numh, listhptr, listh, xijo )
1062
else
1063
nxij = 1
1064
allocate(xijo(3,1))
1065
call memory('A','D',3,'siesta')
1066
endif
1067
C ..................
1068
1069
C Initialize density matrix ...........................................
1070
C set density matrix for first step
1071
found = .false.
1072
dminit = .false.
1073
if (istp .eq. 1) dminit = .true.
1074
if (istp .ne. 1 .and. harrisfun) dminit = .true.
1075
if (istp .ne. 1 .and. (idyn .eq. 6)
1076
$ .and. usesavedm .and. writedm) dminit = .true.
1077
if (istp .ne. 1 .and. (idyn .eq. 7)
1078
$ .and. usesavedm) dminit = .true.
1079
1080
C If auxiliary cell has changed, optionally reset density matrix
1081
C and set usesavedata to false to avoid reading back saved copy
1082
if (initdmaux.and.auxchanged) then
1083
dminit = .true.
1084
usesavedmloc = .false.
1085
else
1086
usesavedmloc = usesavedm
1087
endif
1088
1089
if (dminit)
1090
. call initdm(Datm, Dscf, Dold, lasto, na_s,
1091
. maxnh, no_s, no_l, nspin, na_u, no_l, nspin,
1092
. numh, numhold, listhptr, listhptrold,
1093
. listh, listhold, iaorb, found, inspn,
1094
. usesavedmloc, no_u)
1095
1096
1097
C Initialize energy-density matrix to zero for first call to overfsm
1098
Escf(1:maxnh,1:nspin) = 0.0_dp
1099
1100
C Extrapolate density matrix between steps
1101
itest = .false.
1102
istpsave = 0
1103
iord = 1
1104
if (idyn .eq. 0) iord = 0
1105
if (idyn .eq. 6) iord = 0
1106
if (idyn .eq. 7) iord = 0
1107
C If DM has just been read from disk,
1108
C call extrapol with istep = 2 and iord = 0
1109
C to make it update the structure of DM, if needed
1110
if (found .and. ((istp .eq. 1) .or. (idyn .eq. 6)
1111
. .or. (idyn .eq. 7))) then
1112
istpsave = istp
1113
istp = 2
1114
iord = 0
1115
itest = .true.
1116
endif
1117
if (.not.harrisfun)
1118
. call extrapol(istp, iord, nspin, no_s, no_l, maxnh,
1119
. numh, listhptr, listh, numhold, listhptrold,
1120
. listhold, Dscfsave, Dscf)
1121
C If DM have just been read, restore istp
1122
if (itest) istp = istpsave
1123
itest = .false.
1124
C ..................
1125
1126
C Check for Pulay auxiliary matrices sizes ...................................
1127
if (pulfile .or. maxsav .le. 0) then
1128
nauxpul = 1
1129
if (.not.allocated(auxpul)) then
1130
allocate(auxpul(nauxpul,2))
1131
call memory('A','D',2*nauxpul,'siesta')
1132
endif
1133
else
1134
nauxpul = 0
1135
do io = 1,no_l
1136
nauxpul = nauxpul + numh(io)
1137
enddo
1138
nauxpul = nauxpul * nspin * maxsav
1139
#ifdef MPI
1140
call globalize_max(nauxpul,ntmp)
1141
nauxpul = ntmp
1142
#endif
1143
C Increase nauxpul with safety margin when atoms move
1144
nauxpul = 1.1 * nauxpul + 10
1145
if (allocated(auxpul)) then
1146
if (size(auxpul,1).ne.nauxpul) then
1147
call memory('D','D',size(auxpul),'siesta')
1148
deallocate(auxpul)
1149
allocate(auxpul(nauxpul,2))
1150
call memory('A','D',2*nauxpul,'siesta')
1151
endif
1152
else
1153
allocate(auxpul(nauxpul,2))
1154
call memory('A','D',2*nauxpul,'siesta')
1155
endif
1156
endif
1157
C ....................
1158
1159
C Find overlap matrix ...............................................
1160
call overfsm(na_u, na_s, no_s, scell, xa, indxua, rmaxo, no_l,
1161
. maxna, maxnh, maxnh, lasto, iphorb, isa,
1162
. numh, listhptr, listh, numh, listhptr, listh,
1163
. min(nspin,2), Escf, jna, xij, r2ij,
1164
. fal, stress, S )
1165
C ..................
1166
1167
C Start of SCF iteration _____________________________________________
1168
first = .true.
1169
last = .false.
1170
if (wmix .le. 0._dp) then
1171
if (IOnode) then
1172
write(6,'(/,a,f15.8)')
1173
. 'siesta: WARNING: Mixing weight for SCF loop =', wmix
1174
endif
1175
last = .true.
1176
endif
1177
1178
do iscf = 1, nscf
1179
if (iscf .eq. nscf) last = .true.
1180
call timer( 'IterSCF', 1 )
1181
1182
if (cml_p) call cmlStartStep(xf=mainXML, type='SCF', index=iscf)
1183
1184
C Normalize density matrix to exact charge ...........................
1185
qsol = 0.0_dp
1186
do ispin = 1,min(nspin,2)
1187
do io = 1,nh
1188
qsol = qsol + Dscf(io,ispin) * s(io)
1189
enddo
1190
enddo
1191
#ifdef MPI
1192
call globalize_sum(qsol,buffer1)
1193
qsol = buffer1
1194
#endif
1195
if (IOnode) then
1196
if (.not.first .and.
1197
. abs(qsol/qtot-1._dp).gt.1.d-2) write(6,'(a,2f15.6)')
1198
. 'siesta: WARNING: Qtot, Tr[D*S] =', qtot, qsol
1199
endif
1200
do ispin = 1,nspin
1201
do io = 1,nh
1202
Dscf(io,ispin) = Dscf(io,ispin) * qtot/qsol
1203
Escf(io,ispin) = Escf(io,ispin) * qtot/qsol
1204
enddo
1205
enddo
1206
C ..................
1207
1208
C Initialize Hamiltonian ........................................
1209
H = 0.0_dp
1210
1211
C Initialize forces and stress ...................
1212
if (first.or.last) then
1213
fa(1:3,1:na_u) = 0.0_dp
1214
fal(1:3,1:na_u) = 0.0_dp
1215
stress(1:3,1:3) = 0.0_dp
1216
stressl(1:3,1:3) = 0.0_dp
1217
endif
1218
C ..................
1219
1220
C Self-energy of isolated ions ........................................
1221
if (first) then
1222
Eions = 0.0_dp
1223
do ia = 1,na_u
1224
is = isa(ia)
1225
Eions = Eions + uion(is)
1226
enddo
1227
endif
1228
C ..................
1229
1230
C Neutral-atom: energy, forces and stress ............................
1231
C First time for energy, last time for forces
1232
if (first.or.last) then
1233
call naefs(na_u, na_s, scell, xa, indxua, rmaxv,
1234
. maxna, isa, jna, xij, r2ij,
1235
. Ena, fa, stress)
1236
call dnaefs(na_u, na_s, scell, xa, indxua, rmaxv,
1237
. maxna, isa, jna, xij, r2ij,
1238
. DEna, fa, stress)
1239
Ena = Ena + DEna
1240
endif
1241
C ..................
1242
1243
C Kinetic: energy, forces, stress and matrix elements .................
1244
if (first.or.last) then
1245
call kinefsm(na_u, na_s, no_s, scell, xa, indxua, rmaxo, no_l,
1246
. maxna, maxnh, maxnh, lasto, iphorb, isa,
1247
. numh, listhptr, listh, numh, listhptr, listh,
1248
. min(nspin,2), Dscf, jna, xij, r2ij,
1249
. Ekin, fal, stressl, H )
1250
#ifdef MPI
1251
C Global reduction of energy terms
1252
call globalize_sum(Ekin,buffer1)
1253
Ekin = buffer1
1254
#endif
1255
endif
1256
C ..................
1257
1258
C Non-local-pseudop: energy, forces, stress and matrix elements .......
1259
if (first.or.last) then
1260
call nlefsm(scell, na_u, na_s, isa, xa, indxua, maxna,
1261
. maxnh, maxnh, lasto, lastkb, iphorb, iphKB,
1262
. numh, listhptr, listh, numh, listhptr, listh,
1263
. min(nspin,2), Dscf, Enl, fal, stressl, H)
1264
1265
C Check whether maxna has been increased during nlefsm & resize if needed
1266
if (maxna.gt.size(jna)) then
1267
call memory('D','I',size(jna),'siesta')
1268
deallocate(jna)
1269
call memory('D','D',size(r2ij),'siesta')
1270
deallocate(r2ij)
1271
call memory('D','D',size(xij),'siesta')
1272
deallocate(xij)
1273
allocate(jna(maxna))
1274
call memory('A','I',maxna,'siesta')
1275
allocate(r2ij(maxna))
1276
call memory('A','D',maxna,'siesta')
1277
allocate(xij(3,maxna))
1278
call memory('A','D',3*maxna,'siesta')
1279
endif
1280
1281
#ifdef MPI
1282
C Global reduction of energy terms
1283
call globalize_sum(Enl,buffer1)
1284
Enl = buffer1
1285
#endif
1286
endif
1287
C ..................
1288
1289
C Save or get partial Hamiltonian (non-SCF part) ......................
1290
if (first.or.last) then
1291
do io = 1,nh
1292
H0(io) = H(io,1)
1293
enddo
1294
else
1295
do ispin = 1,nspin
1296
if (ispin .le. 2) then
1297
do io = 1,nh
1298
H(io,ispin) = H0(io)
1299
enddo
1300
else
1301
do io = 1,nh
1302
H(io,ispin) = 0.0_dp
1303
enddo
1304
endif
1305
enddo
1306
endif
1307
C ..................
1308
1309
C Non-SCF part of total energy .......................................
1310
if (first.or.last) then
1311
E0 = -Eions + Ena + Ekin + Enl
1312
else
1313
E0 = 0.0_dp
1314
do ispin = 1,min(nspin,2)
1315
do io = 1,nh
1316
E0 = E0 + H0(io) * Dscf(io,ispin)
1317
enddo
1318
enddo
1319
#ifdef MPI
1320
C Global reduction of E0
1321
call globalize_sum(E0,buffer1)
1322
E0 = buffer1
1323
#endif
1324
E0 = E0 - Eions + Ena
1325
endif
1326
C ..................
1327
1328
C Non-local-pseudop: energy, forces, stress and matrix elements .......
1329
C Add SCF contribution to energy and matrix elements ..................
1330
g2max = g2cut
1331
if (last) then
1332
C Last call to dhscf and grid-cell sampling if requested
1333
ifa = 1
1334
istr = 1
1335
call grdsam( nspin, no_s, iaorb, iphorb,
1336
. no_l, no_u, na_u, na_s, isa, xa, indxua,
1337
. ucell, mscell, g2max, ntm, ifa, istr, maxnh,
1338
. maxnh, numh, listhptr, listh, Dscf, Datm, H,
1339
. Enaatm, Enascf, Uatm, Uscf, DUscf, DUext,
1340
. Exc, Dxc, dipol, fa, stress, fal, stressl)
1341
else
1342
ifa = 0
1343
istr = 0
1344
ihmat = 1
1345
call dhscf( nspin, no_s, iaorb, iphorb, no_l,
1346
. no_u, na_u, na_s, isa, xa, indxua,
1347
. ucell, mscell, g2max, ntm,
1348
. ifa, istr, ihmat, ' ', ' ', ' ', ' ', ' ', ' ',
1349
. maxnh, numh, listhptr, listh, Dscf, Datm,
1350
. maxnh, numh, listhptr, listh, H,
1351
. Enaatm, Enascf, Uatm, Uscf, DUscf, DUext,
1352
. Exc, Dxc, dipol, fa, stress, fal, stressl)
1353
endif
1354
1355
C Output memory use after first call to dhscf
1356
if (istp.eq.1 .and. iscf.eq.1) call printmemory( 6, 0 )
1357
1358
* if (istp.eq.1 .and. iscf.eq.1) write(6,'(/,a,f10.3,a)')
1359
* . 'siesta: dhscf mesh cutoff =', g2max, ' Ry'
1360
1361
C ..................
1362
1363
C Orthonormalization forces ...........................................
1364
if (last) then
1365
call overfsm(na_u, na_s, no_s, scell, xa, indxua,rmaxo,no_l,
1366
. maxna, maxnh, maxnh, lasto, iphorb, isa,
1367
. numh, listhptr, listh, numh, listhptr, listh,
1368
. min(nspin,2), Escf, jna, xij, r2ij,
1369
. fal, stressl, S )
1370
endif
1371
C ..................
1372
1373
C Metadynamics forces
1374
if (lMetaForce.and.(first.or.last)) then
1375
call meta(xa,na_u,ucell,Emeta,fa,stress,last,last)
1376
endif
1377
1378
C Find entropy ........................................................
1379
C Entropy is returned from the call to diagon. To add to the energy
1380
C the entropy computed from the input charge, here it is assigned to the one
1381
C of the former SCF step
1382
1383
Entropy = 0.0_dp
1384
if (isolve .eq. 0) then
1385
if (istp.gt.1 .or. iscf.gt.1) then
1386
Entropy = Entrop
1387
endif
1388
endif
1389
1390
C Save present density matrix ........................................
1391
do is = 1,nspin
1392
do io = 1,nh
1393
Dold(io,is) = Dscf(io,is)
1394
Eold(io,is) = Escf(io,is)
1395
enddo
1396
enddo
1397
1398
C Save Hamiltonian and overlap matrices ............................
1399
if (savehs) then
1400
call iohs( 'write', gamma, no_u, no_s, nspin, indxuo,
1401
$ maxnh, numh, listhptr, listh, H, S, qtot, temp,
1402
$ xijo )
1403
endif
1404
1405
C Solve eigenvalue problem .........................................
1406
if (.not.last) then
1407
if (isolve .eq. 0) then
1408
call diagon(no_s, nspin, nspin, no_l, maxnh, maxnh, no_u,
1409
. numh, listhptr, listh, numh, listhptr, listh,
1410
. H, S, qtot, fixspin, qs, temp, e1, e2,
1411
. gamma, xijo, indxuo, nkpnt, kpoint, kweight,
1412
. eo, qo, Dscf, Escf, ef, efs, Entrop, no_u,
1413
. occtol, iscf, neigwanted)
1414
Ecorrec = 0.0_dp
1415
!
1416
elseif (isolve .eq. 1) then
1417
if (.not. gamma) call die("Cannot do O(N) with k-points.")
1418
call ordern(usesavelwf,ioptlwf,na_s,no_s,no_l,lasto,isa,qa,
1419
. rcoor,rmaxh,ucell,xa,iscf,istp,ncgmax,etol,eta,
1420
. qtot,maxnh,numh,listhptr,listh,H,S,
1421
. chebef,noeta,rcoorcp,beta,pmax,Dscf,Escf,
1422
. Ecorrec,nspin,qs)
1423
Entrop = 0.0_dp
1424
else
1425
call die('siesta: ERROR: wrong solution method')
1426
endif
1427
1428
C Harris-functional energy ............................................
1429
DEharr = 0.0_dp
1430
do ispin = 1,nspin
1431
C const factor takes into account that there are two nondiagonal
1432
C elements in non-collinear spin density matrix, stored as
1433
C ispin=1 => D11; ispin=2 => D22, ispin=3 => Real(D12);
1434
C ispin=4 => Imag(D12)
1435
const = 1._dp
1436
if (ispin .gt. 2) const = 2._dp
1437
do io = 1,nh
1438
DEharr = DEharr + H(io,ispin) * const *
1439
. ( Dscf(io,ispin) - Dold(io,ispin) )
1440
enddo
1441
enddo
1442
#ifdef MPI
1443
C Global reduction of DEharr
1444
call globalize_sum(DEharr,buffer1)
1445
DEharr = buffer1
1446
#endif
1447
C ..................
1448
1449
C Print populations at each SCF step if requested before mixing ......
1450
1451
if (muldeb) then
1452
write (6,"(/a)")
1453
. 'siesta: Mulliken populations before mixing'
1454
call mulliken( mullipop, nspin, na_u, no_u, maxnh,
1455
. numh, listhptr, listh, S, Dscf, isa,
1456
. lasto, iaorb, iphorb )
1457
endif
1458
C ..................
1459
C Mix input and output energy-density and density matrices ............
1460
C Following line for using and saving the density matrix without mix ..
1461
if (wmix.ne.0._dp) then
1462
C Pulay or Broyden mixing
1463
mmix = mix
1464
iiscf = iscf
1465
if (maxsav .le. 0) then
1466
iiscf = 1
1467
if (iscf .ne. 1) mmix = .true.
1468
endif
1469
if (broyden_maxit == 0) then
1470
call pulayx( pulfile, iiscf, mmix, no_l, no_s, maxnh,
1471
. numh, listhptr, nspin, maxsav, wmix, nkick,
1472
. wmixkick, auxpul(1,1), auxpul(1,2), nauxpul,
1473
. Dscf, Dold, dDmax)
1474
else
1475
call broyden_mixing(iscf, mix, no_l, maxnh,
1476
. numh(1:no_l), listhptr(1:no_l), nspin,
1477
$ wmix, nkick, wmixkick, Dscf, Dold, dDmax)
1478
endif
1479
endif
1480
1481
C Ensure that dDmax is the same on all nodes for convergence test/output
1482
#ifdef MPI
1483
call globalize_max(dDmax,buffer1)
1484
dDmax = buffer1
1485
#endif
1486
C ...................
1487
1488
C Print populations at each SCF step, if requested, after mixing ......
1489
1490
if (muldeb) then
1491
write (6,"(/a)")
1492
. 'siesta: Mulliken populations after mixing'
1493
call mulliken( mullipop, nspin, na_u, no_u, maxnh,
1494
. numh, listhptr, listh, S, Dscf, isa,
1495
. lasto, iaorb, iphorb )
1496
endif
1497
C ..................
1498
1499
C Save density matrix on disk, after mixing ...........................
1500
if (writedm) then
1501
if ((idyn .eq. 6) .or. (idyn .eq. 7)) then
1502
if (istp.eq.1)
1503
. call iodm( 'write', maxnh, no_l, nspin,
1504
. numh, listhptr, listh, Dscf, found )
1505
else
1506
call iodm( 'write', maxnh, no_l, nspin,
1507
. numh, listhptr, listh, Dscf, found )
1508
endif
1509
endif
1510
endif !not last
1511
1512
C Add on force field contribution if required
1513
if (first.or.last) then
1514
call twobody(na_u,xa,isa,ucell,Emm,ifa,fa,istr,stress)
1515
endif
1516
1517
C Print energies ......................................................
1518
DEna = Enascf - Enaatm
1519
Etot = E0 + DEna + DUscf + DUext + Exc + Ecorrec + Emad + Emm +
1520
. Emeta
1521
Eharrs = Etot + DEharr
1522
FreeE = Etot - Temp * Entropy
1523
C Recalculating the energy in the last iter (for gridcellsampling)
1524
C but preserving the value of Eharrs1
1525
if (.not.last) Eharrs1 = Eharrs
1526
1527
if (IOnode.and..not.last) then
1528
call siesta_write_energies()
1529
1530
if (harrisfun) then
1531
write(6,"(/a,f14.6,/)") 'siesta: Eharris(eV) = ', Eharrs/eV
1532
if (cml_p) then
1533
call cmlStartPropertyList(mainXML, title='SCF Cycle')
1534
call cmlAddProperty(xf=mainXML, property=Eharrs/eV,
1535
. units="siestaUnits:eV", dictRef="siesta:Eharrs",
1536
. fmt="(f14.7)")
1537
call cmlEndPropertyList(mainXML)
1538
endif
1539
endif
1540
endif
1541
C ...................
1542
1543
! End of one SCF step - flush stdout
1544
if (ionode) then
1545
call pxfflush(6)
1546
call wallclock("-------------- end of scf step")
1547
endif
1548
1549
C If last iteration, exit SCF loop ....................................
1550
if (last) then
1551
do ispin = 1,nspin
1552
do io = 1,nh
1553
Dscf(io,ispin) = Dold(io,ispin)
1554
Escf(io,ispin) = Eold(io,ispin)
1555
enddo
1556
enddo
1557
if (dumpcharge) then
1558
call plcharge( no_s, na_s, no_u, maxnh, maxna, nspin,
1559
. isa, iphorb, indxuo, lasto,
1560
. scell, nsc, xa, rmaxo, datm )
1561
endif
1562
call timer( 'IterSCF', 2 )
1563
if (cml_p) call cmlEndStep(mainXML)
1564
goto 50
1565
endif
1566
C ...................
1567
1568
C If converged, make last iteration to find forces ....................
1569
dEmax = abs(Etot - Elast)
1570
Elast = Etot
1571
if (require_energy_convergence) then
1572
if (dDmax.lt.dDtol.and.dEmax.lt.dEtol) last = .true.
1573
else
1574
if (dDmax.lt.dDtol) last = .true.
1575
endif
1576
C ...................
1577
1578
call timer( 'IterSCF', 2 )
1579
if (istep.eq.inicoor .and. first) call timer( 'IterSCF', 3 )
1580
first = .false.
1581
if (cml_p) call cmlEndStep(mainXML)
1582
enddo
1583
50 continue
1584
C End of SCF iteration_________________________________________________
1585
1586
C Write final Kohn-Sham Energy ........................................
1587
if (cml_p) call cmlStartPropertyList(mainXML,
1588
. title='Final KS Energy')
1589
if (IOnode) then
1590
if ( .not. harrisfun)
1591
. write(6,"(/a,f14.4)") 'siesta: E_KS(eV) = ', Etot/eV
1592
if (cml_p) call cmlAddProperty(xf=mainXML, property=Etot/eV,
1593
. dictref='siesta:E_KS', units='siestaUnits:eV',
1594
. fmt='(f14.6)')
1595
endif
1596
1597
C Substract egg box effect form energy ................................
1598
if (eggbox_block) then
1599
call eggbox('energy',ucell,na_u,isa,ntm,xa,fa,Etot,
1600
. eggbox_block)
1601
if (IOnode)
1602
. write(6,"(/a,f14.4)") 'siesta: E_KS - E_eggbox = ',Etot/eV
1603
if (cml_p) call cmlAddProperty(xf=mainXML, property=Etot/eV,
1604
. dictref='siesta:E_KS_egg', units='siestaUnits:eV',
1605
. fmt='(f14.6)')
1606
endif
1607
if (cml_p) call cmlEndPropertyList(mainXML)
1608
1609
#ifdef MPI
1610
C Global reduction of forces and stresses
1611
allocate(fatmp(3,na_u))
1612
call memory('A','D',3*na_u,'siesta')
1613
call globalize_sum(stressl(1:3,1:3),stresstmp(1:3,1:3))
1614
call globalize_sum(fal(1:3,1:na_u),fatmp(1:3,1:na_u))
1615
stress(1:3,1:3) = stress(1:3,1:3) + stresstmp(1:3,1:3)
1616
fa(1:3,1:na_u) = fa(1:3,1:na_u) + fatmp(1:3,1:na_u)
1617
call memory('D','D',size(fatmp),'siesta')
1618
deallocate(fatmp)
1619
#else
1620
stress(1:3,1:3) = stress(1:3,1:3) + stressl(1:3,1:3)
1621
fa(1:3,1:na_u) = fa(1:3,1:na_u) + fal(1:3,1:na_u)
1622
#endif
1623
1624
C Substract egg box effect from the forces ............................
1625
if (eggbox_block) then
1626
call eggbox('forces',ucell,na_u,isa,ntm,xa,fa,Etot,eggbox_block)
1627
endif
1628
C ...................
1629
! Compute stress without internal molecular pressure
1630
call remove_intramol_pressure(ucell,stress,na_u,xa,fa,mstress)
1631
1632
C Impose constraints to atomic movements by changing forces ...........
1633
if (RemoveIntraMolecularPressure) then
1634
! Consider intramolecular pressure-removal as another
1635
! kind of constraint
1636
call fixed(ucell,mstress,na_u,isa, amass, xa, fa,
1637
$ cstress, cfa, ntcon )
1638
else
1639
call fixed(ucell,stress,na_u,isa, amass, xa, fa,
1640
$ cstress, cfa, ntcon )
1641
endif
1642
1643
C ...................
1644
1645
C Write atomic forces .................................................
1646
fmax = 0.0_dp
1647
cfmax = 0.0_dp
1648
fres = 0.0_dp
1649
do ix = 1,3
1650
ftot(ix) = 0.0_dp
1651
do ia = 1,na_u
1652
ftem = fa(ix,ia)
1653
cftem = cfa(ix,ia)
1654
ftot(ix) = ftot(ix) + ftem
1655
fres = fres + ftem*ftem
1656
fmax = max( fmax, dabs(ftem) )
1657
cfmax = max( cfmax, dabs(cftem) )
1658
enddo
1659
enddo
1660
fres = dsqrt( fres / (3.0_dp*na_u) )
1661
1662
C Calculate and output Zmatrix forces
1663
if (lUseZmatrix) then
1664
call CartesianForce_to_ZmatForce(na_u,xa,fa)
1665
if (IOnode) call iofaZmat
1666
endif
1667
1668
! Compute kinetic contribution to stress
1669
kin_stress(1:3,1:3) = 0.0_dp
1670
do ia = 1,na_u
1671
do jx = 1,3
1672
do ix = 1,3
1673
kin_stress(ix,jx) = kin_stress(ix,jx) -
1674
. amu * amass(ia) * va(ix,ia) * va(jx,ia) / volume
1675
enddo
1676
enddo
1677
enddo
1678
! Add kinetic term to stress tensor
1679
tstress = stress + kin_stress
1680
1681
C Force output .......................................................
1682
if (IOnode) then
1683
call siesta_write_forces()
1684
call siesta_write_stress_pressure()
1685
call wallclock('--- end of geometry step')
1686
endif
1687
1688
C Mulliken population analysis .......................................
1689
call mulliken( mullipop, nspin, na_u, no_u, maxnh,
1690
. numh, listhptr, listh, S, Dscf, isa,
1691
. lasto, iaorb, iphorb )
1692
!
1693
! Save structural information in crystallographic format,
1694
! and canonical Zmatrix (if applicable) before moving the atoms
1695
!
1696
if (IONode) then
1697
call write_struct( ucell, na_u, isa, iza, xa )
1698
if (lUseZmatrix) then
1699
call write_canonical_ucell_and_Zmatrix()
1700
endif
1701
endif
1702
1703
C Save the last coordinates for which the density matrix has been calculated
1704
C or at every coor step if calculating the BECs (the polarisation uses xalast)
1705
if ( (istep.eq.fincoor) .or. bornz ) then
1706
xalast(1:3,1:na_s)=xa(1:3,1:na_s)
1707
endif
1708
1709
Ekinion = 0.0_dp
1710
vn = 0.0_dp
1711
vpr = 0.0_dp
1712
kn = 0.0_dp
1713
kpr = 0.0_dp
1714
1715
iunit = 2
1716
1717
C Move atoms ..........................................................
1718
select case(idyn)
1719
case(0)
1720
if (nmove .ne. 0) then
1721
if (RelaxCellOnly) then
1722
call cell_broyden_optimizer( na_u, xa, ucell,
1723
$ cstress, volume, tp, strtol,
1724
$ varcel, relaxd)
1725
else
1726
1727
if (lUseZmatrix) then
1728
if (broyden_optim) then
1729
call zm_broyden_optimizer( na_u, xa, cfa, ucell,
1730
$ cstress, volume, dxmax, tp, ftol, strtol,
1731
$ varcel, relaxd)
1732
else
1733
call cgvc_zmatrix( na_u, xa, cfa, ucell, cstress,
1734
$ volume, dxmax, tp, ftol, strtol, varcel,
1735
$ relaxd, usesavecg )
1736
endif
1737
else
1738
if (broyden_optim) then
1739
call broyden_optimizer( na_u, xa, cfa, ucell,
1740
$ cstress, volume, dxmax, tp, ftol, strtol,
1741
$ varcel, relaxd )
1742
else
1743
call cgvc( na_u, xa, cfa, ucell, cstress, volume,
1744
$ dxmax, tp, ftol, strtol, varcel,
1745
$ relaxd, usesavecg )
1746
endif
1747
endif
1748
endif ! RelaxCellOnly
1749
1750
! Propagate the new structure to the virtual supercell
1751
call superx( ucell, nsc, na_u, na_s, xa, scell )
1752
if (relaxd) goto 60
1753
! Exit coordinate relaxation loop
1754
endif
1755
!----------
1756
case(1)
1757
call verlet2(istp, iunit, iquench, na_u, cfa, dt,
1758
. amass, ntcon, va, xa, Ekinion, tempion)
1759
! Propagate the new structure to the virtual supercell
1760
call superx( ucell, nsc, na_u, na_s, xa, scell )
1761
! Check convergence for quenching runs (which are
1762
! really relaxations)
1763
if ( iquench .ne. 0 ) then
1764
relaxd = .true.
1765
do ia = 1, na_u
1766
do i = 1, 3
1767
relaxd = relaxd .and. ( abs(cfa(i,ia)) .lt. ftol )
1768
enddo
1769
enddo
1770
endif
1771
!-----------
1772
case (2)
1773
call nose(istp, iunit, na_u, cfa, tt, dt, amass, mn,
1774
. ntcon, va, xa, Ekinion, kn, vn, tempion)
1775
! Propagate the new structure to the virtual supercell
1776
call superx( ucell, nsc, na_u, na_s, xa, scell )
1777
!-----------
1778
case (3)
1779
call pr(istp, iunit, iquench, na_u, cfa, cstress, tp, dt,
1780
. amass, mpr, ntcon, va, xa, vcell, ucell, Ekinion,
1781
. kpr, vpr, tempion, Pint)
1782
! Propagate the new structure to the virtual supercell
1783
call superx( ucell, nsc, na_u, na_s, xa, scell )
1784
if (IOnode) write(6,'(/,a,f12.3,a)')
1785
. 'siesta: E_kin PR =', kpr/Kelvin, ' K'
1786
!-----------
1787
case (4)
1788
call npr(istp, iunit, na_u, cfa, cstress, tp, tt, dt,
1789
. amass, mn, mpr, ntcon, va, xa, vcell, ucell,
1790
. Ekinion, kn, kpr, vn, vpr, tempion, Pint)
1791
! Propagate the new structure to the virtual supercell
1792
call superx( ucell, nsc, na_u, na_s, xa, scell )
1793
!-----------
1794
case (5)
1795
call anneal(istp, iunit, ianneal, taurelax, bulkm,
1796
. na_u, cfa, cstress, tp, tt, dt, amass, ntcon,
1797
. va, xa, ucell, Ekinion, tempion, Pint)
1798
! Propagate the new structure to the virtual supercell
1799
call superx( ucell, nsc, na_u, na_s, xa, scell )
1800
!-----------
1801
case (6:7)
1802
continue !We can't go until after ioxv - see below
1803
!-----------
1804
case (8)
1805
call forcesToPipe( na_u, Etot, cfa, cstress )
1806
end select
1807
1808
if (IOnode) then
1809
if (idyn .gt. 0 .and. idyn .lt. 6) then
1810
write(6,'(/,a,f12.3,a)')
1811
. 'siesta: Temp_ion =', tempion, ' K'
1812
endif
1813
endif
1814
1815
C Save last atomic positions and velocities
1816
C (it should be before moving atoms!)
1817
call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, foundxv)
1818
if (lUseZmatrix)
1819
. call iozm('write',ucell,vcell,xa,foundzm)
1820
call siesta_write_positions
1821
C ...................
1822
1823
C Restore original coordinates after FC displacements
1824
if (idyn .eq. 6 .and. istep .ne. 0) then
1825
xa(ixdispl,iadispl)=xa(ixdispl,iadispl)-dx
1826
endif
1827
if (idyn .eq. 7) then
1828
call phonon_restore_coords(istep,xa,ucell)
1829
endif
1830
1831
C Save atomic positions and velocities accumulatively ................
1832
if (writmd.and.IOnode) then
1833
if ( .not. harrisfun) then
1834
Etot_output = Etot
1835
else
1836
Etot_output = Eharrs1
1837
endif
1838
getot = Etot_output + Ekinion + kn + kpr + vn + vpr
1839
call iomd( na_u, isa, iza,
1840
. xa, va, ucell, vcell, varcel, istep, inicoor,
1841
. fincoor, tempion, Etot_output, getot,
1842
$ volume/Ang**3, Psol/kbar)
1843
call md_v_format(na_u,isa,xa,ucell)
1844
#ifdef CDF
1845
call md_netcdf( na_u, isa, iza,
1846
. xa, va, ucell, vcell, varcel,
1847
. tempion, Etot_output, getot,
1848
$ volume/Ang**3, Psol/kbar)
1849
#endif
1850
1851
endif
1852
1853
C Accumulate coor in Xmol file for animation .........................
1854
lastst = fincoor .le. istep
1855
if (writpx.and.IOnode)
1856
. call pixmol(iza, xa, na_u, lastst)
1857
C ...................
1858
1859
C******************Born charge calculation***************************
1860
if (bornz) then
1861
if (mod(istep,2) .eq. 0 ) then
1862
if (nkpol.lt.1) then
1863
if (IOnode) write(6,'(/,a,f12.6)')
1864
. 'siesta: specify polarization grid for BC calculation'
1865
if (IOnode) write(6,'(a,f12.6)')
1866
. 'siesta: The Born charge matrix will not be calculated'
1867
goto 80
1868
endif
1869
if (IOnode) write(6,'(/,a,f12.6)')
1870
. 'siesta: Calculating polarization. '
1871
1872
C Find total population of spin up and down
1873
if (nspin .ge. 2) then
1874
do ispin = 1,nspin
1875
qspin(ispin) = 0.0_dp
1876
do io = 1,no_l
1877
do j = 1,numh(io)
1878
ind = listhptr(io) + j
1879
qspin(ispin) = qspin(ispin)
1880
. + Dscf(ind,ispin)*S(ind)
1881
enddo
1882
enddo
1883
enddo
1884
#ifdef MPI
1885
C Global reduction of spin components
1886
call globalize_sum(qspin(1:nspin),qtmp(1:nspin))
1887
qspin(1:nspin) = qtmp(1:nspin)
1888
#endif
1889
endif
1890
if (nkpol.gt.0) then
1891
call KSV_pol(na_u, na_s, xalast, rmaxo, scell, ucell,
1892
. no_u, no_l, no_s, nspin, qspin, maxna,
1893
. maxnh, nkpol, numh, listhptr, listh,
1894
. H, S, H0, xijo, indxuo, isa, iphorb,
1895
. iaorb, lasto, jna, xij, r2ij,shape,
1896
. nkpol,kpol,wgthpol, polR, polxyz)
1897
endif
1898
if (nkpol.gt.0.and.IOnode) then
1899
call obc( polxyz, polR, ucell, dx, nspin, node )
1900
endif
1901
endif
1902
endif
1903
80 continue
1904
C*************End born charge calculation******************
1905
1906
60 continue
1907
C Output memory use at the end of this geometry step
1908
if (cml_p) call cmlEndStep(mainXML)
1909
call printmemory( 6, 0 )
1910
call timer( 'IterMD', 2 )
1911
1912
C End of one MD step - flush stdout
1913
if (ionode) call pxfflush(6)
1914
1915
if (relaxd) exit
1916
1917
enddo
1918
C End of coordinate-relaxation loop ==================================
1919
final = .true.
1920
1921
! We want xalast to equal xa for coordinate relaxation only.
1922
if (idyn==0) xalast(1:3,1:na_s)=xa(1:3,1:na_s)
1923
1924
if (cml_p) then
1925
call cmlStartModule(xf=mainXML, title='Finalization')
1926
endif
1927
1928
if (IOnode) then
1929
C Print atomic coordinates (and also unit cell for ParrRah.)
1930
if (nmove .ne. 0) then
1931
if (relaxd)
1932
. call outcoor(ucell, xa, na_u, 'Relaxed', .true. )
1933
if (.not.relaxd)
1934
. call outcoor(ucell, xa, na_u,
1935
. 'Final (unrelaxed)', .true. )
1936
endif
1937
call siesta_write_positions()
1938
if (lUseZmatrix) call write_Zmatrix
1939
if ( varcel .or. (idyn.eq.8)) call outcell(ucell)
1940
1941
C Print coordinates in xmol format in a separate file
1942
1943
if (fdf_boolean('WriteCoorXmol',.false.))
1944
. call coxmol(iza, xa, na_u )
1945
1946
C Print coordinates in cerius format in a separate file
1947
1948
if (fdf_boolean('WriteCoorCerius',.false.))
1949
. call coceri(iza, xa, ucell, na_u, sname )
1950
1951
C Find interatomic distances (output in file BONDS_FINAL)
1952
1953
call bonds( ucell, na_u, isa, xa,
1954
$ rmax_bonds, trim(slabel) // ".BONDS_FINAL" )
1955
1956
endif ! IONode
1957
1958
C Find and print wavefunctions at selected k-points
1959
if (nwk.gt.0) then
1960
call wwave( no_s, nspin, nspin, no_u, no_l, maxnh, maxwk,
1961
. numh, listhptr, listh, H, S, Ef, xijo, indxuo,
1962
. nwk, wfk, no_u, gamma, occtol )
1963
endif
1964
1965
C Find and print bands
1966
if (nbk.gt.0) then
1967
call bands( no_s, nspin, nspin, no_u, no_l, maxnh, maxbk,
1968
. numh, listhptr, listh, H, S, Ef, xijo, indxuo,
1969
. .true., nbk, bk, ebk, no_u, occtol )
1970
if (IOnode) then
1971
if ( writbk ) then
1972
write(6,'(/,a,/,a4,a12)')
1973
. 'siesta: Band k vectors (Bohr**-1):', 'ik', 'k'
1974
do ik = 1,nbk
1975
write(6,'(i4,3f12.6)') ik, (bk(ix,ik),ix=1,3)
1976
enddo
1977
endif
1978
1979
if ( writb ) then
1980
write(6,'(/,a,/,a4,a3,a7)')
1981
. 'siesta: Band energies (eV):', 'ik', 'is', 'eps'
1982
do ispin = 1,min(nspin,2)
1983
do ik = 1,nbk
1984
write(6,'(i4,i3,10f7.2)')
1985
. ik, ispin, (ebk(io,ispin,ik)/eV,io=1,min(10,no_u))
1986
if (no_u.gt.10) write(6,'(7x,10f7.2)')
1987
. (ebk(io,ispin,ik)/eV,io=11,no_u)
1988
enddo
1989
enddo
1990
endif
1991
endif
1992
endif
1993
1994
C Print eigenvalues
1995
if (IOnode .and. writeig) then
1996
if (isolve.eq.0 .and. no_l.lt.1000) then
1997
if (nspin .le. 2) then
1998
write(6,'(/,a,/,a4,a3,a7)')
1999
. 'siesta: Eigenvalues (eV):', 'ik', 'is', 'eps'
2000
do ik = 1,nkpnt
2001
do ispin = 1,nspin
2002
write(6,'(i4,i3,10f7.2)')
2003
. ik,ispin,(eo(io,ispin,ik)/eV,io=1,min(10,neigwanted))
2004
if (no_u.gt.10) write(6,'(7x,10f7.2)')
2005
. (eo(io,ispin,ik)/eV,io=11,neigwanted)
2006
enddo
2007
enddo
2008
else
2009
write(6,'(/,a)') 'siesta: Eigenvalues (eV):'
2010
do ik = 1,nkpnt
2011
write(6,'(a,i6)') 'ik =', ik
2012
write(6,'(10f7.2)')
2013
. ((eo(io,ispin,ik)/eV,io=1,neigwanted),ispin=1,2)
2014
enddo
2015
endif
2016
write(6,'(a,f15.6,a)') 'siesta: Fermi energy =', ef/eV, ' eV'
2017
endif
2018
endif
2019
2020
if (isolve.eq.0.and.IOnode)
2021
. call ioeig(eo,ef,neigwanted,nspin,nkpnt,no_u,nspin,maxk,
2022
. kpoint, kweight)
2023
2024
C Compute the projected density of states
2025
if (IOnode) then
2026
do_pdos = fdf_block('ProjectedDensityOfStates',iu)
2027
if (isolve.ne.0.and.do_pdos) then
2028
write(6,*)
2029
. 'siesta: ERROR: PDOS implemented only with diagon'
2030
do_pdos = .false.
2031
endif
2032
endif
2033
call broadcast(do_pdos)
2034
2035
if (do_pdos) then
2036
C Find the desired energy range
2037
if (IOnode) then
2038
read(iu,'(a)') line
2039
p=>digest(line)
2040
if (nvalues(p).lt.3 .or. nnames(p).ne.1)
2041
$ call die("Wrong format in PDOS block")
2042
factor = fdf_convfac( names(p,1), 'Ry' )
2043
e1 = values(p,1) * factor
2044
e2 = values(p,2) * factor
2045
sigma = values(p,3) * factor
2046
nhist = integers(p,1)
2047
write(6,'(a)') 'siesta: PDOS info: '
2048
write(6,'(a,3(f8.2,a),2x,i5)')
2049
$ 'siesta: e1, e2, sigma, nhist: ',
2050
$ e1/eV,' eV',e2/eV,' eV',sigma/eV,' eV', nhist
2051
endif
2052
2053
call broadcast(e1)
2054
call broadcast(e2)
2055
call broadcast(sigma)
2056
call broadcast(nhist)
2057
2058
call pdos( no_s, nspin, nspin, no_l, maxnh,
2059
. no_u, numh, listhptr, listh, H, S,
2060
. e1, e2, sigma, nhist,
2061
. gamma, xijo, indxuo, nkpnt, kpoint, kweight, eo,
2062
. no_u)
2063
2064
endif ! PDOS calc (do_pdos)
2065
2066
C Print program's energy decomposition and final forces
2067
if (IOnode) then
2068
call siesta_write_energies()
2069
call siesta_write_forces()
2070
call siesta_write_stress_pressure()
2071
call write_basis_enthalpy(FreeE)
2072
endif
2073
2074
C Print spin polarization
2075
if (nspin .ge. 2) then
2076
do ispin = 1,nspin
2077
qspin(ispin) = 0.0_dp
2078
do io = 1,no_l
2079
do j = 1,numh(io)
2080
ind = listhptr(io)+j
2081
jo = listh(ind)
2082
qspin(ispin) = qspin(ispin) + Dscf(ind,ispin) * S(ind)
2083
enddo
2084
enddo
2085
enddo
2086
2087
#ifdef MPI
2088
C Global reduction of spin components
2089
call globalize_sum(qspin(1:nspin),qtmp(1:nspin))
2090
qspin(1:nspin) = qtmp(1:nspin)
2091
#endif
2092
if (nspin .eq. 2) then
2093
if (IOnode) then
2094
write(6,'(/,a,f12.6)')
2095
. 'siesta: Total spin polarization (Qup-Qdown) =',
2096
. qspin(1) - qspin(2)
2097
endif
2098
if (cml_p) call cmlAddProperty(xf=mainXML,
2099
. property=qspin(1)-qspin(2), dictref='siesta:qspin')
2100
elseif (nspin .eq. 4) then
2101
call spnvec( nspin, qspin, qaux, stot, svec )
2102
if (IOnode) then
2103
write(6,'(/,a,f12.6)')
2104
. 'siesta: Total spin polarization (Qup-Qdown) =', stot
2105
write(6,'(a,3f12.6)') 'siesta: Spin vector =', svec
2106
if (cml_p) then
2107
call cmlAddProperty(xf=mainXML, property=stot,
2108
. dictref='siesta:stot')
2109
call cmlAddProperty(xf=mainXML, property=svec,
2110
. dictref='siesta:svec')
2111
endif !cml_p
2112
endif
2113
endif
2114
endif
2115
2116
C Print electric dipole
2117
if (shape .ne. 'bulk') then
2118
if (IOnode) then
2119
write(6,'(/,a,3f12.6)')
2120
. 'siesta: Electric dipole (a.u.) =', dipol
2121
write(6,'(a,3f12.6)')
2122
. 'siesta: Electric dipole (Debye) =',
2123
. (dipol(ix)/Debye,ix=1,3)
2124
endif
2125
if (cml_p) then
2126
call cmlAddProperty(xf=mainXML, property=dipol,
2127
. title='Electric dipole', dictref='siesta:dipol',
2128
. units='siestaUnits:atomic')
2129
endif !cml_p
2130
endif
2131
2132
C Calculation of the bulk polarization using the Berry phase
2133
C formulas by King-Smith and Vanderbilt
2134
C Attention H0 is used as an auxiliary array
2135
if (nkpol.gt.0 .and. .not.bornz) then
2136
call KSV_pol(na_u, na_s, xalast, rmaxo, scell, ucell,
2137
. no_u, no_l, no_s, nspin, qspin, maxna,
2138
. maxnh, nkpol, numh, listhptr, listh,
2139
. H, S, H0, xijo, indxuo, isa, iphorb,
2140
. iaorb, lasto, jna, xij, r2ij,shape,
2141
. nkpol,kpol,wgthpol, polR, polxyz )
2142
endif
2143
2144
C Calculation of the optical conductivity
2145
C Attention H0, Eold, Dold are used as auxiliary arrays
2146
call optical(na_u, na_s, xa, scell, ucell,
2147
. no_u, no_l, no_s, nspin, qspin,
2148
. maxna, maxnh, numh, listhptr, listh, H, S, H0,
2149
. Eold(1,1), Dold(1,1),
2150
. xijo, indxuo, indxua, ebk, ef, temp,
2151
. isa, iphorb, iphKB, iaorb, lasto, lastkb,
2152
. jna, xij, r2ij, shape )
2153
2154
c...................................
2155
2156
C Save electron density and potential
2157
call fdf_global_get(savrho,'SaveRho',
2158
$ dumpcharge .or. .false.)
2159
call fdf_global_get(savdrh,'SaveDeltaRho', .false.)
2160
call fdf_global_get(savevh,'SaveElectrostaticPotential',
2161
$ .false.)
2162
call fdf_global_get(savevt,'SaveTotalPotential', .false.)
2163
call fdf_global_get(savepsch,'SaveIonicCharge', .false.)
2164
call fdf_global_get(savetoch,'SaveTotalCharge', .false.)
2165
2166
if (savrho .or. savdrh .or. savevh .or. savevt .or.
2167
. savepsch .or. savetoch ) then
2168
filrho = ' '
2169
fildrh = ' '
2170
filevh = ' '
2171
filevt = ' '
2172
filepsch = ' '
2173
filetoch = ' '
2174
if (savrho) filrho = paste( slabel, '.RHO' )
2175
if (savdrh) fildrh = paste( slabel, '.DRHO' )
2176
if (savevh) filevh = paste( slabel, '.VH' )
2177
if (savevt) filevt = paste( slabel, '.VT' )
2178
if (savepsch) filepsch = paste( slabel, '.IOCH' )
2179
if (savetoch) filetoch = paste( slabel, '.TOCH' )
2180
g2max = g2cut
2181
call dhscf( nspin, no_s, iaorb, iphorb, no_l,
2182
. no_u, na_u, na_s, isa, xa, indxua,
2183
. ucell, mscell, g2max, ntm,
2184
. 0, 0, 0, filrho, fildrh, filevh, filevt,
2185
. filepsch, filetoch,
2186
. maxnh, numh, listhptr, listh, Dscf, Datm,
2187
. maxnh, numh, listhptr, listh, H,
2188
. Enaatm, Enascf, Uatm, Uscf, DUscf, DUext, Exc, Dxc,
2189
. dipol, fa, stress, fal, stressl )
2190
endif
2191
2192
c Find local density of states
2193
if (IOnode) then
2194
genlogic = fdf_block('LocalDensityOfStates',iu)
2195
endif
2196
call broadcast(genlogic)
2197
2198
if ( genlogic ) then
2199
2200
C Find the desired energy range
2201
if (IOnode) then
2202
read(iu,'(a)') line
2203
p=>digest(line)
2204
if (.not. match(p,"vvn"))
2205
. call die("Wrong format in LocalDensityofStates")
2206
factor = fdf_convfac( names(p,1), 'Ry' )
2207
e1 = values(p,1)*factor
2208
e2 = values(p,2)*factor
2209
call destroy(p)
2210
endif
2211
call broadcast(e1)
2212
call broadcast(e2)
2213
2214
! Find the density matrix for states between e1 and e2
2215
if (isolve .eq. 0) then
2216
call diagon(no_s, nspin, nspin, no_l, maxnh, maxnh, no_u,
2217
. numh, listhptr, listh, numh, listhptr, listh,
2218
. H, S, qtot, fixspin, qs, temp, e1, e2,
2219
. gamma, xijo, indxuo, nkpnt, kpoint, kweight,
2220
. eo, qo, Dscf, Escf, ef, efs, Entrop, no_u,
2221
. occtol, iscf, neigwanted)
2222
2223
! Find the LDOS in the real space mesh
2224
filrho = paste( slabel, '.LDOS' )
2225
g2max = g2cut
2226
call dhscf( nspin, no_s, iaorb, iphorb, no_l,
2227
. no_u, na_u, na_s, isa, xa, indxua,
2228
. ucell, mscell, g2max, ntm,
2229
. 0, 0, 0, filrho, ' ', ' ', ' ', ' ', ' ',
2230
. maxnh, numh, listhptr, listh, Dscf, Datm,
2231
. maxnh, numh, listhptr, listh, H,
2232
. Enaatm, Enascf, Uatm, Uscf, DUscf, DUext, Exc, Dxc,
2233
. dipol, fa, stress, fal, stressl )
2234
else
2235
if (IOnode) write(6,*)
2236
. 'siesta: ERROR: LDOS implemented only with diagon'
2237
endif
2238
2239
endif ! genlogic
2240
2241
C Output memory use up to the end of the program
2242
call printmemory( 6, 1 )
2243
2244
C Print allocation report
2245
call alloc_report( printNow=.true. )
2246
2247
C Stop time counter
2248
call timer( 'siesta', 2 )
2249
call timer( 'all', 3 )
2250
2251
C Print final date and time
2252
if (IOnode) then
2253
call timestamp('End of run')
2254
call wallclock('End of run')
2255
endif
2256
2257
C Finalize MPI
2258
#ifdef MPI
2259
call MPI_Finalize( MPIerror )
2260
#endif
2261
2262
if (cml_p) then
2263
call cmlEndModule(mainXML)
2264
call siesta_cml_exit()
2265
endif
2266
2267
! End of program
2268
! Internal subroutines follow.
2269
2270
contains
2271
2272
2273
subroutine siesta_write_forces()
2274
2275
! Almost the same forces output whether during simulation
2276
! or at the end. Unfortunately not quite, therefore slightly
2277
! tortuous logic below. If we are content to change format
2278
! of output file slightly, this can be simplified.
2279
if (.not.final) then
2280
! print forces to xml every step.
2281
! output forces to stdout depending on writef
2282
if (cml_p) then
2283
call cmlStartPropertyList(mainXML, title='Forces')
2284
call cmlAddProperty(xf=mainXML, property=fa*Ang/eV,
2285
. dictref='siesta:forces', units='siestaUnits:evpa')
2286
call cmlAddProperty(xf=mainXML, property=ftot,
2287
. dictref='siesta:ftot')
2288
call cmlAddProperty(xf=mainXML, property=fmax,
2289
. dictref='siesta:fmax')
2290
call cmlAddProperty(xf=mainXML, property=fres,
2291
. dictref='siesta:fres')
2292
call cmlAddProperty(xf=mainXML, property=cfmax,
2293
. dictref='siesta:cfmax')
2294
call cmlEndPropertyList(mainXML)
2295
endif
2296
write(6,'(/,a)') 'siesta: Atomic forces (eV/Ang):'
2297
if (writef) then
2298
write(6,'(i6,3f12.6)')(ia,(fa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
2299
else
2300
call iofa( na_u, fa )
2301
endif
2302
write(6,'(40("-"),/,a6,3f12.6)') 'Tot',(ftot(ix)*Ang/eV,ix=1,3)
2303
write(6,'(40("-"),/,a6, f12.6)') 'Max',fmax*Ang/eV
2304
write(6,'(a6,f12.6,a)')'Res',fres*Ang/eV,
2305
. ' sqrt( Sum f_i^2 / 3N )'
2306
write(6,'(40("-"),/,a6, f12.6,a)') 'Max',cfmax*Ang/eV,
2307
. ' constrained'
2308
else !not final
2309
C In finalization, only print forces if sufficiently large.
2310
fmax = maxval(abs(fa))
2311
ftot = sum(fa, dim=2)
2312
if (fmax .gt. ftol) then
2313
write(6,'(/,a)') 'siesta: Atomic forces (eV/Ang):'
2314
write(6,'(a,i6,3f12.6)')
2315
. ('siesta: ', ia,(fa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
2316
write(6,'(a,40("-"),/,a,a6,3f12.6)')
2317
. 'siesta: ','siesta: ','Tot',(ftot(ix)*Ang/eV,ix=1,3)
2318
if (cml_p) then
2319
call cmlStartPropertyList(mainXML, title='Force Summary')
2320
call cmlAddProperty(xf=mainXML, property=fa*Ang/eV,
2321
. dictref='siesta:forces', units='siestaUnits:evpa')
2322
call cmlAddProperty(xf=mainXML, property=ftot*Ang/eV,
2323
. dictref='siesta:ftot', units='siestaUnits:evpa')
2324
call cmlEndPropertyList(mainXML)
2325
endif !cml_p
2326
endif
2327
if (Any(cfa /= fa)) then
2328
fmax = maxval(abs(cfa))
2329
ftot = sum(cfa, dim=2)
2330
if (fmax .gt. ftol) then
2331
write(6,'(/,a)') 'siesta: Constrained forces (eV/Ang):'
2332
write(6,'(a,i6,3f12.6)')
2333
. ('siesta: ',ia,(cfa(ix,ia)*Ang/eV,ix=1,3),ia=1,na_u)
2334
write(6,'(a,40("-"),/,a,a4,3f12.6)')
2335
. 'siesta: ','siesta: ','Tot',(ftot(ix)*Ang/eV,ix=1,3)
2336
if (cml_p) then
2337
call cmlStartPropertyList(mainXML,
2338
. title='Constrained Force Summary')
2339
call cmlAddProperty(xf=mainXML, property=cfa*Ang/eV,
2340
. dictref='siesta:cforces', units='siestaUnits:evpa')
2341
call cmlAddProperty(xf=mainXML, property=ftot*Ang/eV,
2342
. dictref='siesta:cftot', units='siestaUnits:evpa')
2343
call cmlEndPropertyList(mainXML)
2344
endif !cml_p
2345
endif
2346
endif
2347
endif !final for forces
2348
2349
end subroutine siesta_write_forces
2350
2351
2352
subroutine siesta_write_stress_pressure()
2353
use zmatrix, only: nZmol, nZmolStartAtom
2354
use zmatrix, only: nZmolAtoms
2355
2356
integer :: ifirst, ilast, natoms_mol, imol
2357
! Stress tensor and pressure:
2358
2359
2360
if (.not.final) then
2361
!
2362
! Write Voigt components of total stress tensor
2363
!
2364
ps = stress + kin_stress
2365
write(6,'(/,a,6f12.2))')
2366
. 'Stress-tensor-Voigt (kbar):',
2367
. (ps(jx,jx)/kbar,jx=1,3),
2368
$ ps(1,2)/kbar,
2369
$ ps(2,3)/kbar,
2370
$ ps(1,3)/kbar
2371
Press = - ((ps(1,1) + ps(2,2) + ps(3,3))/3.0_dp)
2372
write(6,"(a,f14.4)") "Enthalpy (eV/cell)",
2373
$ (FreeE + Press*volume)/eV
2374
2375
if (RemoveIntraMolecularPressure) then
2376
ps = mstress + kin_stress
2377
write(6,'(/,a,6f12.2))')
2378
. 'Non-Molecular-Stress-Voigt (kbar):',
2379
. (ps(jx,jx)/kbar,jx=1,3),
2380
$ ps(1,2)/kbar,
2381
$ ps(2,3)/kbar,
2382
$ ps(1,3)/kbar
2383
Press = - ((ps(1,1) + ps(2,2) + ps(3,3))/3.0_dp)
2384
write(6,"(a,f14.4)") "Non-Molecular Enthalpy (eV/cell)",
2385
$ (FreeE + Press*volume)/eV
2386
endif
2387
!
2388
! Write "target enthalpy" (E + pV, where p is the *target* pressure)
2389
write(6,"(a,f14.4)") "Target enthalpy (eV/cell)",
2390
$ (FreeE + tp*volume)/eV
2391
2392
! Output depends on dynamics option
2393
select case (idyn)
2394
case(0:5,8)
2395
if (idyn==0 .and. (.not.varcel)) then
2396
continue
2397
else
2398
write(6,'(/,a,3(/,a,3f12.6))')
2399
. 'siesta: Stress tensor (static) (eV/Ang**3):',
2400
. (' ',(stress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
2401
Psol = - ((stress(1,1) + stress(2,2) + stress(3,3))/3.0_dp)
2402
write(6,'(/,a,f20.8,a)')
2403
. 'siesta: Pressure (static):', Psol/kBar, ' kBar'
2404
if (cml_p) then
2405
call cmlAddProperty(xf=mainXML, property=stress*Ang**3,
2406
. dictref='siesta:stress')
2407
call cmlAddProperty(xf=mainXML, property=Psol,
2408
. dictref='siesta:psol', title='Pressure (Static)')
2409
endif !cml_p
2410
!
2411
write(6,'(/,a,3(/,a,3f12.6))')
2412
. 'siesta: Stress tensor (total) (eV/Ang**3):',
2413
. (' ',(tstress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
2414
Psol = - ((tstress(1,1)+tstress(2,2) +tstress(3,3))/3.0_dp)
2415
write(6,'(/,a,f20.8,a)')
2416
. 'siesta: Pressure (total):', Psol/kBar, ' kBar'
2417
if (cml_p) then
2418
call cmlAddProperty(xf=mainXML, property=tstress*Ang**3,
2419
. dictref='siesta:tstress')
2420
call cmlAddProperty(xf=mainXML, property=Psol,
2421
. dictref='siesta:tpsol', title='Pressure (Total)')
2422
endif !cml_p
2423
2424
if (RemoveIntraMolecularPressure) then
2425
ps = mstress
2426
write(6,'(/,a,3(/,a,3f12.6))')
2427
. 'siesta: Stress tensor (nonmol) (eV/Ang**3):',
2428
. (' ',(ps(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
2429
Psol = - ((ps(1,1) + ps(2,2) + ps(3,3))/3.0_dp)
2430
write(6,'(/,a,f20.8,a)')
2431
. 'siesta: Pressure (nonmol):', Psol/kBar, ' kBar'
2432
if (cml_p) then
2433
call cmlAddProperty(xf=mainXML, property=ps*Ang**3,
2434
. dictref='siesta:mstress')
2435
call cmlAddProperty(xf=mainXML, property=Psol,
2436
. dictref='siesta:pmol', title='Pressure (Nonmol)')
2437
endif !cml_p
2438
!
2439
ps = mstress + kin_stress
2440
write(6,'(/,a,3(/,a,3f12.6))')
2441
. 'siesta: Stress tensor (nonmol+kin) (eV/Ang**3):',
2442
. (' ',(ps(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
2443
Psol = - ((ps(1,1)+ps(2,2) +ps(3,3))/3.0_dp)
2444
write(6,'(/,a,f20.8,a)')
2445
. 'siesta: Pressure (nonmol+kin):', Psol/kBar, ' kBar'
2446
if (cml_p) then
2447
call cmlAddProperty(xf=mainXML, property=ps*Ang**3,
2448
. dictref='siesta:tmstress')
2449
call cmlAddProperty(xf=mainXML, property=Psol,
2450
. dictref='siesta:tpmol', title='Pressure (Nonmol+Kin)')
2451
endif !cml_p
2452
endif ! Remove intramolecular pressure
2453
2454
endif !varcel
2455
2456
! Write Force Constant matrix if FC calculation ...
2457
case(6)
2458
call ofc(fa,dx,na_u)
2459
case(7)
2460
call phonon_write_forces(fa,na_u,ns,ucell,istep)
2461
end select !idyn
2462
2463
else !final
2464
2465
C Print stress tensor unconditionally
2466
write(6,'(/,a,3(/,a,3f12.6))')
2467
. 'siesta: Stress tensor (static) (eV/Ang**3):',
2468
. ('siesta: ',(stress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
2469
if (cml_p) then
2470
call cmlAddProperty(xf=mainXML, property=stress*Ang**3/eV,
2471
. dictref='siesta:stress', units='siestaUnits:eV_Ang__3')
2472
endif !cml_p
2473
2474
C Print constrained stress tensor if different from unconstrained
2475
if (Any(cstress /= stress )) then
2476
write(6,'(/,a,3(/,a,3f12.6))')
2477
. 'siesta: Constrained stress tensor (static) (eV/Ang**3):',
2478
. ('siesta: ',(cstress(jx,ix)*Ang**3/eV,jx=1,3),ix=1,3)
2479
if (cml_p) then
2480
call cmlAddProperty(xf=mainXML, property=cstress*Ang**3/eV,
2481
. dictref='siesta:cstress',
2482
. units='siestaUnits:eV_Ang__3')
2483
endif !cml_p
2484
endif
2485
2486
C Find pressure
2487
Psol = - (( stress(1,1) + stress(2,2) + stress(3,3) )/3.0_dp)
2488
Pmol = - (( mstress(1,1) + mstress(2,2) + mstress(3,3) )/3.0_dp)
2489
!
2490
write(6,'(/,a,f18.6,a)')
2491
. 'siesta: Cell volume =', volume/Ang**3, ' Ang**3'
2492
write(6,'(/,a,/,a,2a20,a,3(/,a,2f20.8,a))')
2493
. 'siesta: Pressure (static):',
2494
. 'siesta: ','Solid', 'Molecule', ' Units',
2495
. 'siesta: ', Psol, Pmol, ' Ry/Bohr**3',
2496
. 'siesta: ', Psol*Ang**3/eV, Pmol*Ang**3/eV, ' eV/Ang**3',
2497
. 'siesta: ', Psol/kBar, Pmol/kBar, ' kBar'
2498
if (cml_p) then
2499
call cmlStartPropertyList(mainXML, title='Final Pressure')
2500
call cmlAddProperty(xf=mainXML, property=volume/Ang**3,
2501
. title='cell volume', dictref='siesta:cellvol',
2502
. units='siestaUnits:Ang__3')
2503
call cmlAddProperty(xf=mainXML, property=Psol/kBar,
2504
. title='Pressure of Solid', dictref='siesta:pressSol',
2505
. units='siestaUnits:kbar')
2506
call cmlAddProperty(xf=mainXML, property=Pmol/kBar,
2507
. title='Pressure of Molecule', dictref='siesta:pressMol',
2508
. units='siestaUnits:kbar')
2509
call cmlEndPropertyList(mainXML)
2510
endif !cml_p
2511
2512
endif !final for stress & pressure
2513
2514
end subroutine siesta_write_stress_pressure
2515
2516
2517
subroutine siesta_write_energies()
2518
! Only print out full decomposition at very beginning and end.
2519
if ((istp==1.and.first).or.final) then
2520
write(6,'(/,a,/,(a,f17.6))')
2521
. 'siesta: Program''s energy decomposition (eV):',
2522
. 'siesta: Eions =', Eions/eV,
2523
. 'siesta: Ena =', Ena/eV,
2524
. 'siesta: Ekin =', Ekin/eV,
2525
. 'siesta: Enl =', Enl/eV,
2526
. 'siesta: DEna =', DEna/eV,
2527
. 'siesta: DUscf =', DUscf/eV,
2528
. 'siesta: DUext =', DUext/eV,
2529
. 'siesta: Exc =', Exc/eV,
2530
. 'siesta: eta*DQ =', Ecorrec/eV,
2531
. 'siesta: Emadel =', Emad/eV,
2532
. 'siesta: Emeta =', Emeta/eV,
2533
. 'siesta: Emolmec =', Emm/eV,
2534
. 'siesta: Ekinion =', Ekinion/eV,
2535
. 'siesta: Eharris =', (Eharrs1+Ekinion)/eV,
2536
. 'siesta: Etot =', (Etot+Ekinion)/eV,
2537
. 'siesta: FreeEng =', (FreeE+Ekinion)/eV
2538
if (cml_p) then
2539
call cmlStartPropertyList(mainXML,
2540
. title='Energy Decomposition')
2541
call cmlAddProperty(xf=mainXML, property=Eions/eV,
2542
. units='siestaUnits:eV',
2543
. dictref='siesta:Eions', fmt='(f17.6)')
2544
call cmlAddProperty(xf=mainXML, property=Ena/eV,
2545
. units='siestaUnits:eV',
2546
. dictref='siesta:Ena', fmt='(f17.6)')
2547
call cmlAddProperty(xf=mainXML, property=Ekin/eV,
2548
. units='siestaUnits:eV',
2549
. dictref='siesta:Ekin', fmt='(f17.6)')
2550
call cmlAddProperty(xf=mainXML, property=Enl/eV,
2551
. units='siestaUnits:eV',
2552
. dictref='siesta:Enl', fmt='(f17.6)')
2553
call cmlAddProperty(xf=mainXML, property=DEna/eV,
2554
. units='siestaUnits:eV',
2555
. dictref='siesta:DEna', fmt='(f17.6)')
2556
call cmlAddProperty(xf=mainXML, property=DUscf/eV,
2557
. units='siestaUnits:eV',
2558
. dictref='siesta:DUscf', fmt='(f17.6)')
2559
call cmlAddProperty(xf=mainXML, property=DUext/eV,
2560
. units='siestaUnits:eV',
2561
. dictref='siesta:DUext', fmt='(f17.6)')
2562
call cmlAddProperty(xf=mainXML, property=Exc/eV,
2563
. units='siestaUnits:eV',
2564
. dictref='siesta:Exc', fmt='(f17.6)')
2565
call cmlAddProperty(xf=mainXML,property=Ecorrec/eV,
2566
. units='siestaUnits:eV',
2567
. dictref='siesta:Ecorrec', fmt='(f17.6)')
2568
call cmlAddProperty(xf=mainXML, property=Emad/eV,
2569
. units='siestaUnits:eV',
2570
. dictref='siesta:Emad', fmt='(f17.6)')
2571
call cmlAddProperty(xf=mainXML, property=Emeta/eV,
2572
. units='siestaUnits:eV',
2573
. dictref='siesta:Emeta', fmt='(f17.6)')
2574
call cmlAddProperty(xf=mainXML, property=Emm/eV,
2575
. units='siestaUnits:eV',
2576
. dictref='siesta:Emm', fmt='(f17.6)')
2577
call cmlAddProperty(xf=mainXML,property=Ekinion/eV,
2578
. units='siestaUnits:eV',
2579
. dictref='siesta:Ekinion', fmt='(f17.6)')
2580
call cmlAddProperty(xf=mainXML, property=(Eharrs1+Ekinion)/eV,
2581
. units='siestaUnits:eV',
2582
. dictref='siesta:EharrsK', fmt='(f17.6)')
2583
call cmlAddProperty(xf=mainXML, property=(Etot+Ekinion)/eV,
2584
. units='siestaUnits:eV',
2585
. dictref='siesta:EtotK', fmt='(f17.6)')
2586
call cmlAddProperty(xf=mainXML, property=(FreeE+Ekinion)/eV,
2587
. units='siestaUnits:eV',
2588
. dictref='siesta:FreeEK', fmt='(f17.6)')
2589
call cmlEndPropertyList(mainXML)
2590
endif
2591
endif
2592
! On all SCF steps, print out the current energy (format depending on type of run)
2593
if (.not.final) then
2594
! Print total energy and density matrix error .........................
2595
if (cml_p) then
2596
call cmlStartPropertyList(mainXML, title='SCF Cycle')
2597
! Eharrs is always output
2598
call cmlAddProperty(xf=mainXML, property=Eharrs/eV,
2599
. units="siestaUnits:eV",
2600
. dictRef="siesta:Eharrs", fmt="(f14.7)")
2601
endif
2602
! This chain of if statements determines which properties are output.
2603
if (harrisfun) then
2604
write(6,"(/a,f14.6,/)") 'siesta: Eharris(eV) = ', Eharrs/eV
2605
! No need for further cml output
2606
elseif (isolve==0) then
2607
if (cml_p)
2608
. call cmlAddProperty(xf=mainXML, property=FreeE/eV,
2609
. units="siestaUnits:eV",
2610
. dictRef="siesta:FreeE", fmt="(f14.7)")
2611
if (fixspin) then
2612
if (cml_p) then
2613
call cmlAddProperty(xf=mainXML, property=Etot/eV,
2614
. units="siestaUnits:eV",
2615
. dictRef="siesta:Etot", fmt="(f14.7)")
2616
call cmlAddProperty(xf=mainXML, property=FreeE/eV,
2617
. units="siestaUnits:eV",
2618
. dictRef="siesta:FreeE", fmt="(f14.7)")
2619
call cmlAddProperty(xf=mainXML, property=dDmax/eV,
2620
. units="siestaUnits:eV",
2621
. dictRef="siesta:dDmax", fmt="(f14.7)")
2622
endif
2623
if ((iscf .eq. 1).or.muldeb)
2624
. write(6,'(/,a12,3a14,a8,a7,a11)')
2625
. 'siesta: iscf', ' Eharris(eV)',
2626
. ' E_KS(eV)', ' FreeEng(eV)',
2627
. ' dDmax', ' Ef_up', ' Ef_dn(eV)'
2628
write(6,'(a8,i4,3f14.4,f8.4,2f9.4)')
2629
. 'siesta: ',iscf, Eharrs/eV, Etot/eV, FreeE/eV, dDmax,
2630
. (Efs(i)/eV,i=1,2)
2631
if (cml_p) then
2632
call cmlAddProperty(xf=mainXML, property=Efs(1)/eV,
2633
. units="siestaUnits:eV",
2634
. dictRef="siesta:Efs", fmt="(f14.7)")
2635
call cmlAddProperty(xf=mainXML, property=Efs(2)/eV,
2636
. units="siestaUnits:eV",
2637
. dictRef="siesta:Efs", fmt="(f14.7)")
2638
endif
2639
else !fixspin
2640
if ((iscf .eq. 1).or.muldeb)
2641
. write(6,'(/,a12,3a14,2a8)')
2642
. 'siesta: iscf', ' Eharris(eV)',
2643
. ' E_KS(eV)', ' FreeEng(eV)',
2644
. ' dDmax', ' Ef(eV)'
2645
write(6,'(a8,i4,3f14.4,2f8.4)')
2646
. 'siesta: ',iscf, Eharrs/eV, Etot/eV, FreeE/eV,
2647
. dDmax, Ef/eV
2648
if (cml_p) then
2649
call cmlAddProperty(xf=mainXML,property=Ef/eV,
2650
. units="siestaUnits:eV",
2651
. dictRef="siesta:Ef", fmt="(f14.7)")
2652
endif !cml_p
2653
endif !fixspin
2654
elseif (isolve==1) then
2655
write(6,'(/,a15,i4)') 'siesta: iscf = ',iscf
2656
write(6,'(a14,f15.4,a13,f15.4,a10,f7.4/)')
2657
. 'Eharris(eV) = ',Eharrs/eV,
2658
. ' E_KS(eV) = ',Etot/eV,' dDmax = ',dDmax
2659
if (cml_p) then
2660
call cmlAddProperty(xf=mainXML, property=Etot/eV,
2661
. units="siestaUnits:eV",
2662
. dictRef="siesta:Etot", fmt="(f14.7)")
2663
call cmlAddProperty(xf=mainXML, property=dDmax/eV,
2664
. units="siestaUnits:eV",
2665
. dictRef="siesta:dDmax", fmt="(f14.7)")
2666
endif
2667
endif !harrisfun/isolve
2668
2669
if (cml_p) then
2670
call cmlEndPropertyList(mainXML)
2671
endif
2672
2673
else !final
2674
! Print out additional information in finalization.
2675
2676
write(6,'(/,a)') 'siesta: Final energy (eV):'
2677
write(6,'(a,a15,f15.6)')
2678
. 'siesta: ', 'Kinetic =', Ekin/eV,
2679
. 'siesta: ', 'Hartree =', Uscf/eV,
2680
. 'siesta: ', 'Ext. field =', DUext/eV,
2681
. 'siesta: ', 'Exch.-corr. =', Exc/eV,
2682
. 'siesta: ', 'Ion-electron =', (Enascf+Enl+DUscf-Uscf-Uatm)/eV,
2683
. 'siesta: ', 'Ion-ion =', (Ena+Uatm-Enaatm-Eions)/eV,
2684
. 'siesta: ', 'Ekinion =', Ekinion/eV,
2685
. 'siesta: ', 'Total =', (Etot+Ekinion)/eV
2686
if (cml_p) then
2687
call cmlStartPropertyList(xf=mainXML, title='Final Energy')
2688
call cmlAddProperty(xf=mainXML, property=Ekin/eV,
2689
. units='siestaUnits:eV',
2690
. dictref='siesta:Ekin', fmt='(f17.6)')
2691
call cmlAddProperty(xf=mainXML, property=Uscf/eV,
2692
. units='siestaUnits:eV',
2693
. dictref='siesta:Uscf', fmt='(f17.6)')
2694
call cmlAddProperty(xf=mainXML, property=DUext/eV,
2695
. units='siestaUnits:eV',
2696
. dictref='siesta:DUext', fmt='(f17.6)')
2697
call cmlAddProperty(xf=mainXML, property=Exc/eV,
2698
. units='siestaUnits:eV',
2699
. dictref='siesta:Exc', fmt='(f17.6)')
2700
call cmlAddProperty(xf=mainXML,
2701
. property=(Enascf+Enl+DUscf-Uscf-Uatm)/eV,
2702
. units='siestaUnits:eV',
2703
. dictref='siesta:I-e', fmt='(f17.6)')
2704
call cmlAddProperty(xf=mainXML,
2705
. property=(Ena+Uatm-Enaatm-Eions)/eV,
2706
. units='siestaUnits:eV',
2707
. dictref='siesta:I-I', fmt='(f17.6)')
2708
call cmlAddProperty(xf=mainXML, property=Ekinion/eV,
2709
. units='siestaUnits:eV',
2710
. dictref='siesta:Ekinion', fmt='(f17.6)')
2711
call cmlAddProperty(xf=mainXML, property=(Etot+Ekinion)/eV,
2712
. units='siestaUnits:eV',
2713
. dictref='siesta:Etot', fmt='(f17.6)')
2714
call cmlEndPropertyList(mainXML)
2715
endif !cml_p
2716
endif !final
2717
2718
end subroutine siesta_write_energies
2719
2720
subroutine siesta_write_positions
2721
!
2722
! Write out structural information in "crystallography" format,
2723
! and in "canonical" (see Zmatrix module) zmatrix format.
2724
! By now the structure might have changed in response to forces
2725
! and stresses.
2726
2727
call write_struct( ucell, na_u, isa, iza, xa, moved=.true.)
2728
if (lUseZmatrix) then
2729
call write_canonical_ucell_and_Zmatrix
2730
$ (filename="NEXT_ITER.UCELL.ZMATRIX")
2731
endif
2732
2733
if (cml_p) then
2734
call cmlAddMolecule(xf=mainXML, natoms=na_u, elements=elem,
2735
. refs=cisa, coords=xa/Ang, style='x3', fmt='(f12.6)')
2736
call cmlAddLattice(xf=mainXML, cell=ucell/Ang,
2737
. units='siestaUnits:Ang', dictref='siesta:ucell')
2738
endif
2739
end subroutine siesta_write_positions
2740
2741
end program siesta
12
integer :: istep
13
logical :: relaxd
14
15
!------------------------------------------------------------------------- BEGIN
16
call siesta_init()
17
18
! Begin of coordinate relaxation iteration
19
relaxd = .false.
20
istep = inicoor
21
DO WHILE ((istep.le.fincoor) .AND. (.not. relaxd))
22
23
call siesta_forces( istep )
24
25
call siesta_move( istep, relaxd )
26
if (.not. relaxd) then
27
istep = istep + 1
28
endif
29
30
ENDDO
31
! End of coordinate-relaxation loop
32
33
call siesta_analysis( relaxd )
34
35
call siesta_end()
36
!--------------------------------------------------------------------------- END
37
END program siesta
Loggerhead is a web-based interface for Breezy
Version: 2.0.1