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_cell_fire_optim 9! 10! Fire geometry optimization 11! 12 use precision, only : dp 13 use sys, only: die 14 use alloc 15 use fdf, only: fdf_get 16 17 use m_fire 18 use siesta_options, only: dt 19 20 use m_memory, only: memory, mem_stat 21 use parallel, only : ionode 22 23 implicit none 24 25 public :: cell_fire_optimizer 26 private 27 28 CONTAINS 29 subroutine cell_fire_optimizer( na, xa, cell, stress, 30 $ tp, strtol, varcel, relaxd) 31c *************************************************************************** 32c Fire geometry optimization (Cell Only) 33c 34c Written by Alberto Garcia, April 2007 35 36c ******************************** INPUT ************************************ 37c na : number of atoms (will need to keep fractional coords) 38c real*8 tp : Target pressure 39c logical varcel : true if variable cell optimization 40c *************************** INPUT AND OUTPUT ****************************** 41c real*8 cell(3,3) : Matrix of the vectors defining the unit cell 42c input: current step; output: next step 43c cell(ix,ja) is the ix-th component of ja-th vector 44c real*8 xa(3,na) : atomic coordinates 45c real*8 stress(3,3) : Stress tensor components 46c real*8 strtol : Maximum stress tolerance 47c ******************************** OUTPUT *********************************** 48c logical relaxd : True when converged 49c *************************************************************************** 50 51C 52C Modules 53C 54 55 implicit none 56! Subroutine arguments: 57 58 real(dp), intent(in) :: tp, strtol 59 logical, intent(in) :: varcel 60 real(dp), intent(inout) :: stress(3,3), cell(3,3) 61 integer, intent(in) :: na 62 real(dp), intent(inout) :: xa(3,na) 63 logical, intent(out) :: relaxd 64 65c Internal variables and arrays 66 67 real(dp) :: volume, new_volume, trace 68 69 integer i, j, indi 70 71 real(dp) :: celli(3,3) 72 real(dp) :: stress_dif(3,3) 73 74 real(dp), dimension(:), allocatable :: gxa, gfa 75 real(dp), dimension(:), pointer :: deltamax 76 77 78! Saved internal variables: 79 80 logical, save :: frstme = .true. 81 logical, save :: constant_volume 82 real(dp), save :: initial_volume 83 84 85 real(dp), save :: target_stress(3,3), 86 . precon, 87 . strain(3,3), 88 . cellin(3,3) 89 90 type(fire_t), save :: b 91 integer, save :: numel 92 93 logical, save :: fire_debug 94 real(dp), save :: fire_mass 95 real(dp) :: fire_dt, fire_dt_inc, 96 $ fire_dt_dec, fire_alphamax, 97 $ fire_alpha_factor, fire_dtmax 98 integer :: fire_nmin 99 real(dp), parameter :: dstrain_max = 0.1_dp 100 101 real(dp) :: volcel 102 external :: volcel 103c --------------------------------------------------------------------------- 104 105 volume = volcel(cell) 106 107 if ( frstme ) then 108 fire_mass = fdf_get("MD.FIRE.Mass", 1.0_dp) 109 fire_dt = fdf_get("MD.FIRE.TimeStep", dt) 110 fire_dt_inc = fdf_get("MD.FIRE.TimeInc", FIRE_DEF_dt_inc) 111 fire_dt_dec = fdf_get("MD.FIRE.TimeDec", FIRE_DEF_dt_dec) 112 fire_nmin = fdf_get("MD.FIRE.Nmin", FIRE_DEF_nmin) 113 fire_alphamax = fdf_get("MD.FIRE.AlphaMax", FIRE_DEF_alphamax) 114 fire_alpha_factor = fdf_get("MD.FIRE.AlphaFactor", 115 & FIRE_DEF_alpha_factor) 116 fire_dtmax = fdf_get("MD.FIRE.MaxTimeStep", FIRE_DEF_dtmax) 117 fire_dt = fdf_get("MD.FIRE.TimeStep", dt) 118 fire_debug = fdf_get("MD.FIRE.Debug", .false.) 119 120 121 if (varcel ) then 122 numel = 6 123 else 124 call die("no varcel?") 125 endif 126 if (Ionode) then 127 write(6,'(a,i6)') "Cell_Fire_optim: No of elements: ", 128 $ numel 129 endif 130 131 call fire_setup(b, n=numel, dt=fire_dt, 132 $ debug=fire_debug, 133 $ dt_inc=fire_dt_inc, dt_dec=fire_dt_dec, 134 $ alphamax=fire_alphamax, 135 $ alpha_factor=fire_alpha_factor, 136 $ nmin=fire_nmin) 137 138 if ( varcel ) then 139 140C Check if we want a constant-volume simulation 141 constant_volume = fdf_get("MD.ConstantVolume", .false.) 142 143 call get_target_stress(tp,target_stress) 144 if (constant_volume) target_stress = 0.0_dp 145 146 147C Scale factor for strain variables to share magnitude with coordinates 148C ---- (a length in Bohrs typical of bond lengths ..) ------------------ 149 150 precon = fdf_get('MD.PreconditionVariableCell', 151 . 9.4486344d0,'Bohr') 152 153C Initialize absolute strain and save initial cell vectors ------------- 154C Initialization to 1. for numerical reasons, later substracted -------- 155 156 strain(1:3,1:3) = 1.0_dp 157 cellin(1:3,1:3) = cell(1:3,1:3) 158 initial_volume = volcel(cellin) 159 160 endif 161 162 frstme = .false. 163 endif ! First call 164 165C Variable cell ------------------------------------------------------------- 166 167 if ( varcel ) then 168 169 allocate(gfa(numel), stat=mem_stat) 170 call memory('A','D',numel,'Fire_optim', 171 $ stat=mem_stat,id="gfa") 172 allocate(gxa(numel), stat=mem_stat) 173 call memory('A','D',numel,'Fire_optim', 174 $ stat=mem_stat,id="gxa") 175 nullify( deltamax ) 176 call re_alloc( deltamax, 1, numel, name='deltamax', 177 & routine='fire_optimizer' ) 178 179 180 relaxd = .true. 181 182C Symmetrizing the stress tensor 183 184 do i = 1, 3 185 do j = i+1, 3 186 stress(i,j) = 0.5_dp*( stress(i,j) + stress(j,i) ) 187 stress(j,i) = stress(i,j) 188 enddo 189 enddo 190 191C Subtract target stress 192 193 stress_dif = stress - target_stress 194! 195! Take 1/3 of the trace out here if constant-volume needed 196! 197 if (constant_volume) then 198 trace = stress_dif(1,1) + stress_dif(2,2) + stress_dif(3,3) 199 do i=1,3 200 stress_dif(i,i) = stress_dif(i,i) - trace/3.0_dp 201 enddo 202 endif 203 204C Append excess stress and strain to gxa and gfa ------ 205C preconditioning: scale stress and strain so as to behave similarly to x,f - 206 207 indi = 0 208 do i = 1, 3 209 do j = i, 3 210 indi = indi + 1 211 gfa(indi) = -stress_dif(i,j)*volume/precon 212 gxa(indi) = strain(i,j) * precon 213 deltamax(indi) = dstrain_max 214 enddo 215 enddo 216 217C Check stress convergence -------------------------------------------------- 218 219 do i = 1, 3 220 do j = 1, 3 221 relaxd = relaxd .and. 222 . ( abs(stress_dif(i,j)) .lt. abs(strtol) ) 223 enddo 224 enddo 225 226 if (relaxd) RETURN 227 228C Call Fire step 229 230 call fire_step(b,gfa,gxa,deltamax) 231 232 endif 233 234C Transform back if variable cell 235 236 if ( varcel ) then 237 238 ! New cell 239 240 indi = 0 241 do i = 1, 3 242 do j = i, 3 243 indi = indi + 1 244 strain(i,j) = gxa(indi) / precon 245 strain(j,i) = strain(i,j) 246 enddo 247 enddo 248 249 ! Inverse of matrix of cell vectors (transpose of) 250 call reclat( cell, celli, 0 ) 251 252! Update cell 253 254 cell = cellin + matmul(strain-1.0_dp,cellin) 255 if (constant_volume) then 256 new_volume = volcel(cell) 257 cell = cell * (initial_volume/new_volume)**(1.0_dp/3.0_dp) 258 endif 259 260 call rescale_coordinates(na,xa, celli_oldcell=celli, 261 $ new_cell=cell) 262 263 264 ! Deallocate local memory 265 266 deallocate (gxa, stat=mem_stat) 267 call memory('D','D',numel,'Fire_optim', 268 $ stat=mem_stat,id="gxa") 269 deallocate (gfa, stat=mem_stat) 270 call memory('D','D',numel,'Fire_optim', 271 $ stat=mem_stat,id="gfa") 272 call de_alloc( deltamax, name='deltamax' ) 273 274 endif ! variable cell 275 276 end subroutine cell_fire_optimizer 277! 278!-------------------------------------------------------------- 279 280 subroutine rescale_coordinates(na,xa, 281 $ celli_oldcell,new_cell) 282 use precision, only : dp 283 use zmatrix 284 285 integer, intent(in) :: na 286 real(dp), intent(inout) :: xa(3,na) 287 real(dp), intent(in) :: celli_oldcell(3,3) 288 real(dp), intent(in) :: new_cell(3,3) 289 290 291 real(dp), dimension(3) :: r, frac 292 integer :: ifirst, imol, icart, i, j 293 294 ! NOTE: We have to be careful here if using a Zmatrix 295 if (lUseZmatrix) then 296 297 ! re-scale only the positions of the leading atoms 298 ! in each molecule, 299 ! plus any cartesian block, 300 ! and recompute the cartesian coordinate array 301 ! 302 do imol = 1, nZmol 303 ifirst = nZmolStartAtom(imol) 304 r(1:3) = Zmat(3*(ifirst-1)+1:3*(ifirst-1)+3) 305 frac(1:3) = matmul(transpose(celli_oldcell),r(1:3)) 306 r(1:3) = matmul(new_cell,frac(1:3)) 307 Zmat(3*(ifirst-1)+1:3*(ifirst-1)+3) = r(1:3) 308 enddo 309 do icart = 1, nZcart 310 ! Process cartesian block 311 ifirst = nZcartStartAtom(icart) 312 do j = ifirst, ifirst + nZcartAtoms(icart) - 1 313 r(1:3) = Zmat(3*(j-1)+1:3*(j-1)+3) 314 frac(1:3) = matmul(transpose(celli_oldcell),r(1:3)) 315 Zmat(3*(j-1)+1:3*(j-1)+3) = matmul(new_cell,frac(1:3)) 316 enddo 317 enddo 318 call Zmat_to_Cartesian(xa) 319 else 320 ! No Zmatrix 321 ! Rescale coordinates for all atoms 322 do i = 1, na 323 xa(:,i) = matmul(transpose(celli_oldcell),xa(:,i)) 324 xa(:,i) = matmul(new_cell,xa(:,i)) 325 enddo 326 327 endif 328 329 end subroutine rescale_coordinates 330 331 332 end module m_cell_fire_optim 333 334