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 subroutine cgvc( na, xa, fa, cell, stress, dxmax, 9 . tp, ftol, strtol, varcel, relaxd, usesavecg ) 10c *************************************************************************** 11c Variable-cell conjugate-gradient geometry optimization 12c 13c Energy minimization including atomic coordinates and unit cell vectors. 14c It allows an external target stress: 15c %block MD.TargetStress 16c 3.5 0.0 0.0 0.0 0.0 0.0 17c %endblock MD.TargetStress 18c corresponding to xx, yy, zz, xy, xz, yz. 19c In units of (-MD.TargetPressure) 20c Default: hydrostatic pressure: -1, -1, -1, 0, 0, 0 21c 22c Based on E({xa},stress), with {xa} in fractional coor 23c The gradient of the energy given by {cfa} forces (fractional!) 24c The fractional coordinates are multiplied by the initial cell vectors 25c to get them in Bohr for dxmax and preconditioning. 26c 27c Written by E. Artacho. November 1998. 28c ******************************** INPUT ************************************ 29c integer na : Number of atoms in the simulation cell 30c real*8 fa(3,na) : Atomic forces 31c real*8 dxmax : Maximum atomic (or lattice vector) displacement 32c real*8 tp : Target pressure 33c real*8 ftol : Maximum force tolerance 34c logical varcel : true if variable cell optimization 35c logical usesavecg : true if we're using saved CG files. 36c *************************** INPUT AND OUTPUT ****************************** 37c real*8 xa(3,na) : Atomic coordinates 38c input: current step; output: next step 39c real*8 cell(3,3) : Matrix of the vectors defining the unit cell 40c input: current step; output: next step 41c cell(ix,ja) is the ix-th component of ja-th vector 42c real*8 stress(3,3) : Stress tensor components 43c real*8 strtol : Maximum stress tolerance 44c ******************************** OUTPUT *********************************** 45c logical relaxd : True when converged 46c *************************************************************************** 47 48C 49C Modules 50C 51 use precision, only : dp 52 use parallel, only : Node 53 use sys, only : die 54 use fdf 55 56 use m_mpi_utils, only: broadcast 57 use units, only: kBar, Ang 58! 59! Use old-style conjugate-gradient routine 60! Use the new one only for the Zmatrix case 61! 62 use m_conjgr_old, only: conjgr 63 use alloc, only: re_alloc, de_alloc 64 65 implicit none 66 67! Subroutine arguments: 68 69 integer, intent(in) :: na 70 real(dp), intent(in) :: fa(3,na), dxmax, 71 . tp, ftol, strtol 72 logical, intent(in) :: varcel, usesavecg 73 real(dp), intent(inout) :: xa(3,na), stress(3,3), cell(3,3) 74 logical, intent(out) :: relaxd 75 76c Internal variables and arrays 77 78 real(dp) :: new_volume, trace, ftol_tmp, volume 79 80 logical found 81 integer ia, i, j, n, indi 82 83 real(dp) :: celli(3,3), sxx, syy, szz, sxy, sxz, syz 84 real(dp) :: stress_dif(3,3) 85 86 real(dp), pointer :: gxa(:)=>null(), gfa(:)=>null() 87 real(dp), pointer, save :: cgaux(:)=>null() 88 89 type(block_fdf) :: bfdf 90 type(parsed_line), pointer :: pline=>null() 91 92! Saved internal variables: 93 94 logical, save :: frstme = .true. 95 logical, save :: tarstr = .false. 96 logical, save :: constant_volume 97 real(dp), save :: initial_volume 98 99 real(dp), save :: cgcntr(0:20) = 0.0_dp 100 101 integer, save :: ndeg, 102 . linmin 103 104 real(dp), save :: tstres(3,3), 105 . modcel(3), 106 . precon, 107 . strain(3,3), 108 . cellin(3,3) 109 110 real(dp) :: volcel 111 external :: volcel 112c --------------------------------------------------------------------------- 113 114 volume = volcel(cell) 115 116C Allocate local memory 117 call re_alloc( gfa, 1, na*3+6, 'gfa', 'cgvc' ) 118 call re_alloc( gxa, 1, na*3+6, 'gxa', 'cgvc' ) 119 if (.not.associated(cgaux)) 120 & call re_alloc( cgaux, 1, na*6+12, 'cgaux', 'cgvc' ) 121 122C If first call to cgvc, check dim and get target stress -------------------- 123 124 if ( frstme ) then 125 126 if ( varcel ) then 127 128C Check if we want a constant-volume simulation 129 constant_volume = fdf_get("MD.ConstantVolume", .false.) 130 131C Look for target stress and read it if found, otherwise generate it -------- 132 133 tarstr = fdf_block('MD.TargetStress',bfdf) 134 135 if (tarstr) then 136 if (Node.eq.0) then 137 write(6,'(/a,a)') 'cgvc: Reading %block MD.TargetStress', 138 . ' (units of MD.TargetPressure).' 139 endif 140 if (.not. fdf_bline(bfdf,pline)) 141 . call die('cgvc: ERROR in MD.TargetStress block') 142 sxx = fdf_bvalues(pline,1) 143 syy = fdf_bvalues(pline,2) 144 szz = fdf_bvalues(pline,3) 145 sxy = fdf_bvalues(pline,4) 146 sxz = fdf_bvalues(pline,5) 147 syz = fdf_bvalues(pline,6) 148 call fdf_bclose(bfdf) 149 150 tstres(1,1) = - sxx * tp 151 tstres(2,2) = - syy * tp 152 tstres(3,3) = - szz * tp 153 tstres(1,2) = - sxy * tp 154 tstres(2,1) = - sxy * tp 155 tstres(1,3) = - sxz * tp 156 tstres(3,1) = - sxz * tp 157 tstres(2,3) = - syz * tp 158 tstres(3,2) = - syz * tp 159 else 160 if (Node.eq.0) then 161 write(6,'(/a,a)') 'cgvc: No target stress found, ', 162 . 'assuming hydrostatic MD.TargetPressure.' 163 endif 164 do i= 1, 3 165 do j= 1, 3 166 tstres(i,j) = 0._dp 167 enddo 168 tstres(i,i) = - tp 169 enddo 170 endif 171 172C Write target stress down -------------------------------------------------- 173 174 if (constant_volume) then 175 tstres(:,:) = 0.0_dp 176 if (Node.eq.0) then 177 write(6,"(a)") "***Target stress set to zero " // 178 $ "for constant-volume calculation" 179 endif 180 endif 181 if (Node.eq.0) then 182 write(6,"(/a)") 'cgvc: Target stress (kBar)' 183 write(6,"(a,2x,3f12.3)") 184 . 'cgvc:', tstres(1,1)/kBar, tstres(1,2)/kBar, 185 . tstres(1,3)/kBar 186 write(6,"(a,2x,3f12.3)") 187 . 'cgvc:', tstres(2,1)/kBar, tstres(2,2)/kBar, 188 . tstres(2,3)/kBar 189 write(6,"(a,2x,3f12.3)") 190 . 'cgvc:', tstres(3,1)/kBar, tstres(3,2)/kBar, 191 . tstres(3,3)/kBar 192 endif 193 194 195C Moduli of original cell vectors for fractional coor scaling back to au --- 196 197 do n = 1, 3 198 modcel(n) = 0.0_dp 199 do j = 1, 3 200 modcel(n) = modcel(n) + cell(j,n)*cell(j,n) 201 enddo 202 modcel(n) = dsqrt( modcel(n) ) 203 enddo 204 205C Scale factor for strain variables to share magnitude with coordinates 206C ---- (a length in Bohrs typical of bond lengths ..) ------------------ 207 208 precon = fdf_get('MD.PreconditionVariableCell', 209 . 9.4486344d0,'Bohr') 210 211C Dimension of space where E is minimized ------------------------------ 212 213 ndeg = na*3 + 6 214 215C Initialize absolute strain and save initial cell vectors ------------- 216C Initialization to 1. for numerical reasons, later substracted -------- 217 218 strain(1:3,1:3) = 1.0_dp 219 cellin(1:3,1:3) = cell(1:3,1:3) 220 initial_volume = volcel(cellin) 221 222 else 223 224 ndeg = na*3 225 226 endif 227 228C Initialize and read cgaux and cgcntr if present and wanted --------------- 229 230 if (usesavecg) then 231 if (Node.eq.0) then 232 call iocg( 'read', ndeg*2, cgaux, cgcntr, relaxd, found ) 233 endif 234 call broadcast(cgaux) 235 call broadcast(cgcntr) 236 call broadcast(relaxd) 237 call broadcast(found) 238 239 if ( found ) then 240 linmin = cgcntr(1) 241! Simple bugfix for fixed cell restarts. In the case of a tightend ftol 242! after a sucessful CG run, relaxed is set to .true. by the read from 243! iocg. This results in the test for convergence being skipped within 244! conjgr. Setting relaxed to .false. here avoids this. In the case of 245! a variable cell relaxed is set below anyway and the converence test 246! is always performed. A better solution (moving the test out into a 247! seperate routine) will follow in a later version. -- AMW 8/7/2008 248 relaxd = .false. 249 else 250 if (Node.eq.0) then 251 write(6,'(/,a)') 'cgvc: WARNING: CG file not found' 252 endif 253 relaxd = .false. 254 cgcntr(0) = 0 255 linmin = 1 256 endif 257 else 258 relaxd = .false. 259 cgcntr(0) = 0 260 linmin = 1 261 endif 262 263 frstme = .false. 264 endif 265 266C Variable cell ------------------------------------------------------------- 267 268 if ( varcel ) then 269 270C Inverse of matrix of cell vectors (transpose of) ------------------------ 271 272 call reclat( cell, celli, 0 ) 273 274C Transform coordinates and forces to fractional ---------------------------- 275C but scale them again to Bohr by using the (fix) moduli of the original ---- 276C lattice vectors (allows using maximum displacement as before) ------------- 277C convergence is checked here for input forces as compared with ftol -------- 278 279 relaxd = .true. 280 do ia = 1, na 281 do n = 1, 3 282 indi = 3*(ia - 1) + n 283 gxa(indi) = 0.0_dp 284 gfa(indi) = 0.0_dp 285 relaxd = relaxd .and. ( abs(fa(n,ia)) .lt. ftol ) 286 do i = 1, 3 287 gxa(indi) = gxa(indi) + celli(i,n) * xa(i,ia) * modcel(n) 288 gfa(indi) = gfa(indi) + cell(i,n) * fa(i,ia) / modcel(n) 289 enddo 290 enddo 291 enddo 292 293C Symmetrizing the stress tensor -------------------------------------------- 294 295 do i = 1, 3 296 do j = i+1, 3 297 stress(i,j) = 0.5_dp*( stress(i,j) + stress(j,i) ) 298 stress(j,i) = stress(i,j) 299 enddo 300 enddo 301 302C Subtract target stress 303 304 stress_dif = stress - tstres 305! 306! Take 1/3 of the trace out here if constant-volume needed 307! 308 if (constant_volume) then 309 trace = stress_dif(1,1) + stress_dif(2,2) + stress_dif(3,3) 310 do i=1,3 311 stress_dif(i,i) = stress_dif(i,i) - trace/3.0_dp 312 enddo 313 endif 314 315C Append excess stress and strain to gxa and gfa ------ 316C preconditioning: scale stress and strain so as to behave similarly to x,f - 317 318 indi = 3*na 319 do i = 1, 3 320 do j = i, 3 321 indi = indi + 1 322 gfa(indi) = -stress_dif(i,j)*volume/precon 323 gxa(indi) = strain(i,j) * precon 324 enddo 325 enddo 326 327C Check stress convergence -------------------------------------------------- 328 329 do i = 1, 3 330 do j = 1, 3 331 relaxd = relaxd .and. 332 . ( abs(stress_dif(i,j)) .lt. abs(strtol) ) 333 enddo 334 enddo 335 336C Call conjugate gradient minimization -------------------------------------- 337 338 if ( .not. relaxd ) then 339 ftol_tmp = 0.0_dp 340!! do i=1,ndeg 341!! print *, "gxa, gfa ", i, gxa(i), gfa(i) 342!! enddo 343 call conjgr(ndeg,gxa,gfa,dxmax, ftol_tmp ,cgcntr,cgaux ) 344 endif 345 346 347C Fixed cell ---------------------------------------------------------------- 348 349 else 350 351 if (.not. relaxd) then 352!! do i=1,ndeg 353!! print *, "gxa, gfa ", i, gxa(i), gfa(i) 354!! enddo 355 call conjgr( 3*na, xa, fa, dxmax, ftol, cgcntr, cgaux ) 356 relaxd = (int(cgcntr(0)) .eq. 0) !! ??? 357 endif 358 359 endif 360 361C Checking line minimizations and convergence ------------------------------- 362 363 if (nint(cgcntr(1)) .ne. linmin) then 364 if (Node.eq.0) then 365 write(6,'(/a,i4,a,f10.4)') 366 . 'cgvc: Finished line minimization ', linmin, 367 . '. Mean atomic displacement =', cgcntr(18)/sqrt(dble(na)) 368 endif 369 linmin = nint(cgcntr(1)) 370 endif 371 372C Transform back if variable cell ------------------------------------------- 373 374 if ( varcel .and. (.not. relaxd) ) then 375 376C New cell ------------------------------------------------------------------ 377 378 indi = 3*na 379 do i = 1, 3 380 do j = i, 3 381 indi = indi + 1 382 strain(i,j) = gxa(indi) / precon 383 strain(j,i) = strain(i,j) 384 enddo 385 enddo 386 387 cell = cellin + matmul(strain-1.0_dp,cellin) 388 if (constant_volume) then 389 new_volume = volcel(cell) 390 if (Node.eq.0) write(6,"(a,f12.4)") 391 $ "Volume before coercion: ", new_volume/Ang**3 392 cell = cell * (initial_volume/new_volume)**(1.0_dp/3.0_dp) 393 endif 394 395C Output fractional coordinates to cartesian Bohr, and copy to xa ----------- 396 397 do ia = 1, na 398 do i = 1, 3 399 xa(i,ia) = 0.0_dp 400 do n = 1, 3 401 indi = 3*(ia - 1) + n 402 xa(i,ia) = xa(i,ia) + cell(i,n) * gxa(indi) / modcel(n) 403 enddo 404 enddo 405 enddo 406 407 endif 408 409C Save cgaux ---------------------------------------------------------------- 410 411 if (Node.eq.0) 412 . call iocg( 'write', ndeg*2, cgaux, cgcntr, relaxd, found ) 413 414C Deallocate local memory 415 call de_alloc( gxa, 'gxa', 'cgvc' ) 416 call de_alloc( gfa, 'gfa', 'cgvc' ) 417 418 return 419 end 420 421