1! ---
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt .
6! See Docs/Contributors.txt for a list of contributors.
7! ---
8      module m_struct_init
9      private
10      public :: struct_init
11
12      CONTAINS
13
14      subroutine struct_init()
15      USE siesta_options
16      use units, only: Ang
17      use m_ioxv
18      USE alloc,          only: re_alloc, de_alloc
19
20      use siesta_geom
21      use atmfuncs,       only : izofis
22      use atomlist,       only : elem, iza
23      use fdf
24      use m_mpi_utils,    only : broadcast
25      use periodic_table, only : symbol
26      use m_iostruct,     only : read_struct
27      use siesta_cml
28      use files,          only : slabel
29      use zmatrix,        only : lUseZmatrix
30      use parallel,       only : IOnode
31#ifdef NCDF_4
32      use m_steps, only : inicoor, fincoor
33      use m_exp_coord, only : read_exp_coord_options
34      use m_exp_coord, only : exp_coord_init
35      use m_exp_coord, only : exp_coord_next
36#endif
37      use siesta_master,  only: coordsFromMaster  ! Get coords from master prog
38      use siesta_master,  only: setPropertyForMaster  ! Set a property to be
39                                                      ! sent to master program
40      use siesta_master,  only: siesta_server         ! Is siesta a server?
41
42      implicit none
43
44      real(dp), external :: volcel
45      external :: automatic_cell, bonds, iozm, shaper
46
47      character(len=80) :: dummy_str
48      logical           :: foundxv, foundzm
49      real(dp)          :: local_charnet
50
51      real(dp):: bcell(3,3), tmp_cell(3,3)
52      integer :: i, ia, ix
53      integer:: nbcell  ! Auxiliary to call shaper routine
54      real(dp), pointer :: xang(:,:)
55      character(len=22) :: dyntyp
56      character(len=60) :: restart_file
57      logical :: found_restart, step_back
58
59!------------------------------------------------------------------------- BEGIN
60#ifdef DEBUG
61      call write_debug( '  PRE struct_init' )
62#endif
63
64      !! Read number of atoms and coordinates, and unit cell
65      use_struct_file = fdf_get("UseStructFile",.false.)
66      use_struct_file = fdf_get("MD.UseStructFile",
67     $     use_struct_file)    ! For legacy symbol
68      if (use_struct_file) then
69         call read_struct( na_u, tmp_cell) ! Sets na_u, xa, and isa
70      else
71         call coor(na_u,tmp_cell)  ! Sets na_u, xa, and isa
72      endif
73      ucell = tmp_cell  ! initialize module variable
74
75!     Prepare iza here: it might be needed by ioxv
76      nullify( iza, va )
77      call re_alloc( iza, 1, na_u, 'iza', 'struct_init' )
78      call re_alloc( va, 1, 3, 1, na_u, 'va', 'struct_init' )
79
80      ! Options read here instead of in siesta_options
81      usesavexv = fdf_get('MD.UseSaveXV', usesaveddata)
82      usesavezm = fdf_get('MD.UseSaveZM', usesaveddata)
83      writic    = fdf_get('WriteCoorInitial', .true.)
84
85      ! Read Z-matrix coordinates and forces from file
86      if (lUseZmatrix) then
87        foundzm = .false.
88        if (usesavezm) then
89            call iozm('read',ucell,vcell,xa,foundzm)
90            if (IOnode) then
91                if (.not.foundzm) then
92                   write(6,'(/,a)') 'siesta: WARNING: ZM file not found'
93                else
94                    write(6,'(/,a)')
95     .        "! Info in XV file prevails over previous structure input"
96                endif
97            endif
98        endif
99      endif
100
101! Read cell shape and atomic positions from a former run ..............
102#ifdef NCDF_4
103      dummy_str = fdf_get('MD.TypeOfRun','none')
104      if ( leqi(dummy_str,'explicit') ) then
105         ! Read in file name
106         call read_exp_coord_options( )
107         call exp_coord_init(slabel,na_u,inicoor,fincoor)
108         if ( inicoor > 0 ) then
109            call exp_coord_next(inicoor,na_u,xa)
110         end if
111      end if
112#endif
113      foundxv = .false.
114      if (usesavexv) then
115        call ioxv('read', ucell, vcell, na_u, isa, iza, xa, va, foundxv)
116        if (IOnode) then
117           if (.not.foundxv) then
118              write(6,'(/,a)') 'siesta: WARNING: XV file not found'
119           else
120              write(6,"(a)")
121     $       "! Info in XV file prevails over previous structure input"
122           endif
123        endif
124!       For a Verlet/Nose/PR/NPR run with more than one time step, if
125!       the RESTART file is not found, backward-propagate the atomic
126!       positions by one time step using the Euler method
127
128        istart = fdf_get('MD.InitialTimeStep',1)
129        ifinal = fdf_get('MD.FinalTimeStep',1)
130        if (foundxv .and. (ifinal-istart>0)) then
131           dyntyp = fdf_get('MD.TypeOfRun','CG')
132           step_back=.true.
133           if (leqi(dyntyp,'verlet')) then
134             restart_file = trim(slabel) // '.VERLET_RESTART'
135           else if (leqi(dyntyp,'nose')) then
136             restart_file = trim(slabel) // '.NOSE_RESTART'
137           else if (leqi(dyntyp,'parrinellorahman')) then
138             restart_file = trim(slabel) // '.PR_RESTART'
139           else if (leqi(dyntyp,'noseparrinellorahman')) then
140             restart_file = trim(slabel) // '.NPR_RESTART'
141           else
142             step_back=.false.
143           endif
144           if (step_back) then
145              if (IOnode) then
146                 inquire( file=restart_file, exist=found_restart )
147              endif
148              call broadcast(found_restart)
149              if (.not. found_restart) then
150                 ! dt_default = 1.0 (fs)
151                 ! The user should not rely on the default here...
152                 ! The parameter dt_default is private to read_options
153                 dt = fdf_get('MD.LengthTimeStep',1.0_dp,'fs')
154                 xa=xa-dt*va
155                 if (IOnode) then
156                    write(6,'(5a)') 'WARNING: ', trim(restart_file),
157     $                           ' not found--reading only XV file',
158     $                           ' and moving back 1 time step using',
159     $                           ' Euler'
160                 endif
161              endif
162           endif
163        endif
164      endif
165! ..................
166
167! Find if siesta is a server of forces, unless we already know it .....
168      if (.not.siesta_server) then
169        dummy_str = fdf_get('MD.TypeOfRun',"none")
170        siesta_server = (leqi(dummy_str,'forces') .or.
171     .                    leqi(dummy_str,'master'))
172      end if
173! .....................
174
175! Read cell shape and atomic positions from driver program
176      if (siesta_server) then
177        call setPropertyForMaster('atomic_numbers',
178     .                             na_u, real(iza(:),dp), ' ')
179        call coordsFromMaster( na_u, xa, ucell )
180      end if
181! .....................
182
183! Dump initial coordinates to output ..................................
184      if ( writic.and.(IOnode) ) then
185        write(6,'(/a)') 'siesta: Atomic coordinates (Bohr) and species'
186        write(6,"('siesta: ',2x,3f10.5,i3,3x,i6)")
187     .           ( (xa(ix,ia), ix=1,3), isa(ia), ia, ia=1, na_u)
188      endif
189! ..................
190
191! Automatic cell generation ...........................................
192      if (volcel(ucell) .lt. 1.0d-8) then
193         local_charnet = fdf_get('NetCharge',0.0_dp)
194         call automatic_cell(ucell,scell,na_u,xa,isa,local_charnet)
195      endif
196
197! Initialize atomic velocities to zero ................................
198      if (.not. foundxv) then
199        ! AG ** What happens with iozm call?
200        va(1:3,1:na_u) = 0.0_dp
201        vcell(1:3,1:3) = 0.0_dp
202      endif
203! ..................
204
205
206! Find system shape ...................................................
207      call shaper( ucell, na_u, isa, xa, shape, nbcell, bcell )
208      if (IOnode) then
209        write(6,'(/,2a)') 'siesta: System type = ', shape
210      endif
211
212! Find interatomic distances (output in file BONDS)
213      if (IOnode) then
214        rmax_bonds = fdf_physical("MaxBondDistance", 6.0_dp, "Bohr")
215        call bonds( ucell, na_u, isa, xa,
216     $       rmax_bonds, trim(slabel)// ".BONDS" )
217      endif
218
219! Output of initial system details:
220
221      if (cml_p) then
222! We need the names of the elements on node 0
223        nullify(elem)
224        allocate( elem(na_u) )
225!        call re_alloc( elem, 1, na_u, 'elem', 'struct_init' )
226        do i = 1, na_u
227           elem(i) = symbol(izofis(isa(i)))
228        enddo
229
230        call cmlStartModule(xf=mainXML, title='Initial System')
231        nullify(xang)
232        call re_alloc( xang, 1, 3, 1, na_u, 'xang', 'struct_init' )
233        xang = xa(1:3,1:na_u)/Ang
234        call cmlAddMolecule(xf=mainXML, natoms=na_u,
235     .       coords=xa/Ang, elements=elem, atomRefs=cisa)
236        call de_alloc( xang, 'xang', 'struct_init' )
237
238        call cmlAddLattice(xf=mainXML, cell=ucell,
239     .       units='siestaUnits:angstrom', dictref='siesta:ucell')
240        call cmlAddProperty(xf=mainXML, value=trim(shape),
241     .       dictref='siesta:shape')
242        call cmlEndModule(xf=mainXML)
243      endif
244
245#ifdef DEBUG
246      call write_debug( '  POS struct_init' )
247#endif
248      END subroutine struct_init
249      end module m_struct_init
250