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_siesta_move 9 10 private 11 public :: siesta_move 12 13 CONTAINS 14 15 subroutine siesta_move( istep, relaxd ) 16 17 use zmatrix, only: lUseZmatrix, iofaZmat, write_Zmatrix, 18 . CartesianForce_to_ZmatForce 19 use atomlist, only: iza, amass 20 use m_ioxv, only: ioxv 21 USE siesta_options 22 use units, only: Kelvin 23 use parallel, only: IOnode 24 use siesta_cml, only: cml_p, cmlEndStep, mainXML 25 use sys, only: die 26 27 use m_broyden_optim, only: broyden_optimizer 28 use m_fire_optim, only: fire_optimizer 29 use m_zm_broyden_optim, only: zm_broyden_optimizer 30 use m_zm_fire_optim, only: zm_fire_optimizer 31 use m_cell_broyden_optim, only: cell_broyden_optimizer 32 use m_cell_fire_optim, only: cell_fire_optimizer 33 34 use m_dynamics, only: nose, verlet2, npr, anneal, pr 35 36 use m_energies, only: Ekinion 37 use write_subs 38 use m_steps 39 use m_kinetic 40 use siesta_geom, only: na_u ! Number of atoms in unit cell 41 use siesta_geom, only: na_s ! Number of atoms in super-cell 42 use siesta_geom, only: isa ! Atomic species 43 use siesta_geom, only: xa ! Atomic positions 44 use siesta_geom, only: xa_last ! Previous atomic positions 45 use siesta_geom, only: va ! Atomic velocities 46 use siesta_geom, only: ucell ! Unit cell vectors 47 use siesta_geom, only: ucell_last ! Previous unit cell vectors 48 use siesta_geom, only: scell, nsc ! Super-cell vectors and multip. 49 use siesta_geom, only: scell_last ! Previous super-cell vectors 50 use siesta_geom, only: vcell ! Time derivative of unit cell vecs 51 use m_forces, only: cfa ! Constrained forces 52 use m_forces, only: ntcon ! Number of geometry constraints 53 use m_stress, only: cstress ! Constrained stress tensor 54 use siesta_master, only: coordsFromMaster ! Get coords from master prog 55 use m_target_stress, only: subtract_target_stress 56 57 use m_check_walltime 58 use m_mpi_utils, only: barrier 59 60#ifdef NCDF_4 61 use m_exp_coord, only : exp_coord_next 62#endif 63#ifdef SIESTA__FLOOK 64 use siesta_dicts, only : dict_variable_add 65 use flook_siesta, only : slua_call, LUA_MOVE 66#endif 67 68 implicit none 69 70 integer, intent(in) :: istep 71 logical, intent(out) :: relaxd 72 73 74 integer :: ix, iadispl, ixdispl 75 integer, parameter :: iunit = 2 ! Physical-units option for MD: 76 ! 1=>(eV,Ang), 2=>(Ry,Bohr) 77 real(dp) :: Pint ! Instantaneous pressure 78 real(dp) :: eff_stress(3,3) 79 real(dp) :: tp_pr 80 81 logical :: foundxv ! dummy for call to ioxv 82 logical :: foundzm ! dummy for call to iozm 83 84 logical :: time_is_up 85 character(len=40) :: tmp_str 86 87 real(dp), external :: volcel 88 89#ifdef DEBUG 90 call write_debug( ' PRE siesta_move' ) 91#endif 92 93!------------------------------------------------------------------ BEGIN 94 95 call timer( 'siesta_move', 1 ) 96 97 ! Save the last geometry for which the density matrix 98 ! (electronic structure) has been calculated 99 100 xa_last(1:3,1:na_s) = xa(1:3,1:na_s) 101 ucell_last(1:3,1:3) = ucell(1:3,1:3) 102 scell_last(1:3,1:3) = scell(1:3,1:3) 103 104! Move atoms .......................................................... 105! 106! ** Regarding the output of coordinates: 107! 108! -- 'state_init' should record the actual coordinates of a step (CG, MD, whatever) 109! -- this routine (siesta_move) should only record checkpointing information, such as 110! the predicted coordinates for a possible next geometry step. 111! Historically, the checkpointing has been done with the (highly confusing) 112! XV file (routine ioxv), but there is the newer option (useful only for 113! relaxations) of writing a '.STRUCT_NEXT_ITER' file (routine write_positions with move="true") 114! -- 'siesta_analysis' should only write final coordinates. 115 116! There is a pletora of output routines: (and more in 'siesta_analysis') 117! 118! -- ioxv 119! -- iozm 120! -- write_positions (which calls write_struct and zm_canonical...) 121! -- write_md_record 122! 123! They should be rationalized. 124! 125! In the following, it is better to put all the logic in the individual blocks, avoiding 126! branching and re-checks of the 'idyn' variable (wich, by the way, should have symbolic 127! values instead of hardwired numerical ones) 128 129 select case(idyn) 130 131 case(0) 132 ! The original "relaxation" case, but note that it also covers 133 ! pure single-point calculations (with MD.NumCGsteps = 0) 134 135 136 ! Pure single-point calculations will not enter this block 137 if (nmove .ne. 0) then ! That is, if requesting "CG" steps 138 139 ! Here we want checkpointing, except if the structure is already relaxed 140 141 ! Note that these routines do not update the coordinates if 142 ! the force/stress criterion for relaxation is satisfied 143 144 if (RelaxCellOnly) then 145 if (broyden_optim) then 146 call cell_broyden_optimizer( na_u, xa, ucell, 147 $ cstress, tp, strtol, 148 $ varcel, relaxd) 149 elseif (fire_optim) then 150 call cell_fire_optimizer( na_u, xa, ucell, 151 $ cstress, tp, strtol, 152 $ varcel, relaxd) 153 else 154 call die("Cell-only optim needs Broyden or FIRE") 155 endif 156 157 else ! Coordinate relaxation (and maybe cell) 158 159 if (lUseZmatrix) then 160 if (broyden_optim) then 161 call zm_broyden_optimizer( na_u, xa, ucell, cstress, tp, 162 & strtol, varcel, relaxd ) 163 elseif (fire_optim) then 164 call zm_fire_optimizer( na_u, xa, cfa, ucell, 165 $ cstress, dxmax, tp, ftol, strtol, 166 $ varcel, relaxd) 167 else 168 call cgvc_zmatrix( na_u, xa, cfa, ucell, cstress, 169 $ dxmax, tp, ftol, strtol, varcel, 170 $ relaxd, usesavecg ) 171 endif 172 else 173 if (broyden_optim) then 174 call broyden_optimizer( na_u, xa, cfa, ucell, 175 $ cstress, tp, ftol, strtol, varcel, relaxd ) 176 elseif (fire_optim) then 177 call fire_optimizer( na_u, xa, cfa, ucell, 178 $ cstress, dxmax, tp, ftol, strtol, 179 $ varcel, relaxd ) 180 else 181 call cgvc( na_u, xa, cfa, ucell, cstress, 182 $ dxmax, tp, ftol, strtol, varcel, 183 $ relaxd, usesavecg ) 184 endif 185 endif 186 endif ! RelaxCellOnly 187 188 if (relaxd) then 189 190 ! Will not call ioxv et al in the block below, so that: 191 ! - The XV file will contain the xa coords 192 ! written in the previous step (if any, otherwise 193 ! those written by state_init) 194 ! - The .STRUCT_NEXT_ITER file would be that 195 ! written in the previous step (if any), and it 196 ! will contain the same coords as the .STRUCT_OUT 197 ! file produced in state_init in this step. 198 ! - The Zmatrix files... 199 else 200 call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, 201 & foundxv) 202 if (lUseZmatrix) call iozm('write',ucell,vcell,xa,foundzm) 203 ! This writes the "next_iter" STRUCT and canonical zmatrix files 204 call siesta_write_positions(moved=.true.) 205 ! Save atomic positions and velocities accumulatively and 206 ! accumulate coor in Xmol file for animation 207 call write_md_record( istep ) 208 call superx_if_compat( ucell, nsc, na_u, na_s, xa, scell ) 209 endif 210 211 endif 212 213 case(1) ! Micro-canonical MD with fixed cell (Verlet) 214 215 relaxd = .false. ! It is intent(out) 216 ! Check convergence for quenching runs (which are 217 ! really relaxations) 218 ! Avoid calling the routine if relaxed, so that the coordinates 219 ! are not changed. 220 if ( iquench .ne. 0 ) then 221 relaxd = all(abs(cfa(1:3,1:na_u)) < ftol) 222 ! We might want to fall back on a different relaxation 223 ! scheme when we reach the slow-moving part 224 ! (check temp_ion?; check cfa behavior?) 225 endif 226 if (.not. relaxd) then 227 call verlet2(istp, iunit, iquench, na_u, cfa, dt, 228 . amass, ntcon, va, xa, Ekinion, tempion) 229 if (IOnode) then 230 write(6,'(/,a,f12.3,a)') 231 $ 'siesta: Temp_ion =', tempion, ' K' 232 endif 233 call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, 234 & foundxv) 235 call siesta_write_positions(moved=.true.) 236 call write_md_record( istep ) 237 call superx_if_compat( ucell, nsc, na_u, na_s, xa, scell ) 238 endif 239 240 241 case (2) ! Canonical MD with fixed cell (Nose-Hoover) 242 243 relaxd = .false. ! It is intent(out) 244 call nose( istp, iunit, na_u, cfa, tt, dt, amass, mn, 245 . ntcon, va, xa, Ekinion, kn, vn, tempion ) 246 if (IOnode) then 247 write(6,'(/,a,f12.3,a)') 'siesta: Temp_ion =', tempion, ' K' 248 endif 249 call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, 250 & foundxv) 251 call siesta_write_positions(moved=.true.) 252 call write_md_record( istep ) 253 call superx_if_compat( ucell, nsc, na_u, na_s, xa, scell ) 254 255 256 case (3) ! Micro-canonical variable-cell MD (Parrinello-Rahman) 257 258 ! Check convergence for quenching runs (which are 259 ! really relaxations) 260 ! Avoid calling the routine if relaxed, so that the coordinates 261 ! are not changed. 262 if ( iquench .ne. 0 ) then 263 ! Allow a general target stress in this relaxation mode 264 call subtract_target_stress(cstress,eff_stress) 265 relaxd = all(abs(cfa(1:3,1:na_u)) < ftol) 266 relaxd = relaxd .AND. (all(abs(eff_stress) < strtol)) 267 cstress = eff_stress 268 tp_pr = 0.0_dp ! As we are already passing the modified stress 269 else 270 ! Standard MD mode 271 relaxd = .false. 272 tp_pr = tp ! Keep the original target pressure 273 endif 274 if (.not. relaxd) then 275 call pr(istp, iunit, iquench, na_u, cfa, cstress, tp_pr, dt, 276 . amass, mpr, ntcon, va, xa, vcell, ucell, Ekinion, 277 . kpr, vpr, tempion, Pint) 278 if (IOnode) then 279 write(6,'(/,a,f12.3,a)') 280 . 'siesta: E_kin PR =', kpr/Kelvin, ' K' 281 write(6,'(/,a,f12.3,a)') 282 $ 'siesta: Temp_ion =', tempion, ' K' 283 endif 284 call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, 285 & foundxv) 286 call siesta_write_positions(moved=.true.) 287 call write_md_record( istep ) 288 call superx_if_compat( ucell, nsc, na_u, na_s, xa, scell ) 289 endif 290 291 case (4) ! Canonical variable-cell MD (Nose-Parrinello-Rahman) 292 293 relaxd = .false. 294 call npr(istp, iunit, na_u, cfa, cstress, tp, tt, dt, 295 . amass, mn, mpr, ntcon, va, xa, vcell, ucell, 296 . Ekinion, kn, kpr, vn, vpr, tempion, Pint) 297 if (IOnode) then 298 write(6,'(/,a,f12.3,a)') 299 $ 'siesta: Temp_ion =', tempion, ' K' 300 endif 301 call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, 302 & foundxv) 303 call siesta_write_positions(moved=.true.) 304 call write_md_record( istep ) 305 call superx_if_compat( ucell, nsc, na_u, na_s, xa, scell ) 306 307 case (5) ! Annealings 308 309 relaxd = .false. ! It is intent(out) 310 call anneal(istp, iunit, ianneal, taurelax, bulkm, 311 . na_u, cfa, cstress, tp, tt, dt, amass, ntcon, 312 . va, xa, ucell, Ekinion, tempion, Pint) 313 if (IOnode) then 314 write(6,'(/,a,f12.3,a)') 315 $ 'siesta: Temp_ion =', tempion, ' K' 316 endif 317 call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, 318 & foundxv) 319 call siesta_write_positions(moved=.true.) 320 call write_md_record( istep ) 321 call superx_if_compat( ucell, nsc, na_u, na_s, xa, scell ) 322 323 case (6) ! Force-constant-matrix calculation 324 325 relaxd = .false. ! It is intent(out) 326 327 ! Output of coordinates is meaningless here in general, 328 ! except maybe to checkpoint a calculation to be restarted 329 ! later, by reading the XV file and choosing inicoor and 330 ! fincoor appropriately 331 332 ! Save current (displaced) atomic positions 333 call ioxv( 'write', ucell, vcell, na_u, isa, iza, 334 & xa, va, foundxv) 335 336 ! Undo the last atom displacement 337 if ( istep > inicoor ) then ! inicoor == 0 338 iadispl = (istep-mod(istep-1,6))/6+ia1 339 ix = mod(istep-1,6)+1 340 ixdispl = (ix - mod(ix-1,2) +1)/2 341 xa(ixdispl,iadispl) = xa(ixdispl,iadispl) - dx 342 end if 343 344 ! Displace atom by dx 345 ! NOTE: MOVE is before istep increment, hence 346 ! fincoor == istep will be the last step and 347 ! the initial geometry should be retained. 348 if ( istep < fincoor ) then ! fincoor = (ia2-ia1+1)*6 349 iadispl = ((istep+1)-mod(istep,6))/6+ia1 350 ix = mod(istep,6)+1 351 ixdispl = (ix - mod(ix-1,2) +1)/2 352 dx = -dx 353 xa(ixdispl,iadispl) = xa(ixdispl,iadispl) + dx 354 end if 355 356 case (7) ! PHONON interface --- removed 357 358 relaxd = .true. 359 360 case (8) ! Server mode 361 362 relaxd = .false. 363 364 ! It can be argued that in server operation the 365 ! responsibility of writing coordinates falls on the 366 ! client. 367 368 ! Save atomic positions and velocities accumulatively and 369 ! accumulate coor in Xmol file for animation 370 call write_md_record( istep ) 371 372 ! Get coordinates from driver program 373 call coordsFromMaster( na_u, xa, ucell ) 374 if (volcel(ucell) < 1.0e-8_dp) then 375 call automatic_cell(ucell,scell,na_u,xa,isa,charnet) 376 endif 377 378#ifdef NCDF_4 379 case (9) ! explicit coordinate mode 380 381 relaxd = .false. 382 383 ! Save the md-record 384 call write_md_record( istep ) 385 if ( istep < fincoor ) then 386 call exp_coord_next( istep + 1 , na_u, xa) 387 else 388 relaxd = .true. 389 end if 390#endif 391 392#ifdef SIESTA__FLOOK 393 case (10) ! LUA hook coordinates 394 395 relaxd = .false. 396 397 ! Add the quest for relaxed list of variables 398 call dict_variable_add('MD.Relaxed',relaxd) 399 400 ! Communicate with lua 401 call slua_call(LUA, LUA_MOVE) 402 403 if ( volcel(ucell) < 1.0e-8_dp ) then 404 ! In case the user requests a recalculation of the 405 ! unit-cell size (for molecules) 406 call automatic_cell(ucell,scell,na_u,xa,isa,charnet) 407 end if 408 409 if ( relaxd ) then 410 ! see comments for idyn == 1 411 else 412 call ioxv( 'write', ucell, vcell, na_u, isa, iza, xa, va, 413 & foundxv) 414 415 call siesta_write_positions( moved=.true. ) 416 call write_md_record( istep ) 417 call superx_if_compat( ucell, nsc, na_u, na_s, xa, scell ) 418 end if 419#endif 420 421 end select 422 423! Output memory use at the end of this geometry step 424 if (cml_p) then 425 call cmlEndStep(mainXML) 426 endif 427 428 call timer( 'siesta_move', 2 ) 429 call timer( 'IterGeom', 2 ) 430! End of one MD step - flush stdout 431 if (ionode) call pxfflush(6) 432 433 ! Check whether we are short of time to continue 434 call check_walltime(time_is_up) 435 if (time_is_up) then 436 ! Do any other bookeeping not done by "die" 437 call timer('all',2) 438 call timer('all',3) 439 if (.not. relaxd) then 440 call message('WARNING','GEOM_NOT_CONV: '// 441 $ 'Geometry relaxation not converged') 442 write(tmp_str,"(a,i5)") 'Geometry step: ', istep 443 call message(' (info)', trim(tmp_str)) 444 endif 445 call barrier() 446 call die("OUT_OF_TIME: Time is up at end of geometry step. ") 447 endif 448 449!--------------------------------------------------------------------------- END 450 451#ifdef DEBUG 452 call write_debug( ' POS siesta_move' ) 453#endif 454 455 contains 456 457 ! For backwards compatibility only 458 459 subroutine superx_if_compat( ucell, nsc, na, maxa, xa, scell ) 460 use siesta_options, only : compat_pre_v4_dynamics 461 use atomlist, only : superx 462 463 integer :: maxa, na, nsc(3) 464 real(dp) :: ucell(3,3), xa(3,MAXA), scell(3,3) 465 ! Only xa and scell are overwritten upon exit 466 ! but superx does not have intent (so we ca not either) 467 468 if (compat_pre_v4_dynamics) then 469 call superx( UCELL, NSC, NA, MAXA, XA, SCELL ) 470 end if 471 end subroutine superx_if_compat 472 473 end subroutine siesta_move 474 475 end module m_siesta_move 476