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