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_zm_broyden_optim 9 public :: zm_broyden_optimizer 10 private 11 12 CONTAINS 13! 14 subroutine zm_broyden_optimizer( na, xa, cell, stress, tp, strtol, 15 & varcel, relaxd ) 16c ******************************** INPUT ************************************ 17c integer na : Number of atoms in the simulation cell 18c real*8 stress(3,3) : Stress tensor components 19c real*8 tp : Target pressure 20c real*8 strtol : Maximum stress tolerance 21c logical varcel : true if variable cell optimization 22c *************************** INPUT AND OUTPUT ****************************** 23c real*8 xa(3,na) : Atomic coordinates 24c input: current step; output: next step 25c real*8 cell(3,3) : Matrix of the vectors defining the unit cell 26c input: current step; output: next step 27c cell(ix,ja) is the ix-th component of ja-th vector 28c real*8 stress(3,3) : Stress tensor components 29c real*8 strtol : Maximum stress tolerance 30c ******************************** OUTPUT *********************************** 31c logical relaxd : True when converged 32c *************************************************************************** 33 34C 35C Modules 36C 37 use precision, only : dp 38 use parallel, only : Node, IONode 39 use sys, only : die 40 use fdf 41 use units, only: kBar, Ang 42 use m_broyddj_nocomm, only: broyden_t 43 use m_broyddj_nocomm, only: broyden_reset, broyden_step, 44 $ broyden_destroy, broyden_init, 45 $ broyden_is_setup 46 47 use zmatrix, only : VaryZmat, Zmat 48 use zmatrix, only : CartesianForce_to_ZmatForce 49 use zmatrix, only : ZmatForce,ZmatForceVar 50 use zmatrix, only : iZmattoVars,ZmatType 51 use zmatrix, only : Zmat_to_Cartesian 52 use zmatrix, only : coeffA, coeffB, iNextDept 53 use Zmatrix, only : ZmatForceTolLen, ZmatForceTolAng 54 use Zmatrix, only : ZmatMaxDisplLen, ZmatMaxDisplAng 55 56 implicit none 57 58C Subroutine arguments: 59 integer, intent(in) :: na 60 logical, intent(in) :: varcel 61 logical, intent(out) :: relaxd 62 real(dp), intent(in) :: tp 63 real(dp), intent(inout):: xa(3,na), stress(3,3), strtol, cell(3,3) 64 65C Internal variables and arrays 66 integer :: ia, i, j, n, indi,indi1,vi,k 67 real(dp) :: volume, new_volume, trace 68 real(dp) :: sxx, syy, szz, sxy, sxz, syz 69 real(dp) :: celli(3,3) 70 real(dp) :: stress_dif(3,3) 71 real(dp), allocatable :: gxa(:), gfa(:) 72 real(dp) :: force, force1 73 74C Saved internal variables: 75 integer, save :: ndeg 76 logical, save :: frstme = .true. 77 logical, save :: tarstr = .false. 78 logical, save :: constant_volume 79 real(dp), save :: initial_volume 80 81 real(dp), save :: tstres(3,3) 82 real(dp), save :: cellin(3,3) 83 real(dp), save :: modcel(3) 84 real(dp), save :: precon 85 real(dp), save :: strain(3,3) ! Special treatment !! 86 87 real(dp), allocatable :: ftoln(:), dxmaxn(:), rnew(:) 88 89 type(broyden_t), save :: br 90 logical, save :: initialization_done = .false. 91 92 type(block_fdf) :: bfdf 93 type(parsed_line), pointer :: pline 94 95 real(dp), save :: jinv0 96 integer, save :: maxit 97 logical, save :: cycle_on_maxit, variable_weight 98 logical, save :: broyden_debug 99 real(dp) :: weight 100 101 integer :: ndegi, ndi 102 103 real(dp), external :: volcel 104 105c --------------------------------------------------------------------------- 106 107 volume = volcel(cell) 108 109 if (.not. initialization_done) then 110 111 maxit = fdf_get("MD.Broyden.History.Steps",5) 112 cycle_on_maxit = fdf_get("MD.Broyden.Cycle.On.Maxit",.true.) 113 variable_weight = fdf_get("MD.Broyden.Variable.Weight",.false.) 114 broyden_debug = fdf_get("MD.Broyden.Debug",.false.) 115 jinv0 = fdf_get("MD.Broyden.Initial.Inverse.Jacobian", 116 $ 1.0_dp) 117 118 if (ionode) then 119 write(6,*) 120 write(6,'(a,i3)') 121 $ "Broyden_optim: max_history for broyden: ", maxit 122 write(6,'(a,l1)') 123 $ "Broyden_optim: cycle on maxit: ", cycle_on_maxit 124 if (variable_weight) write(6,'(a)') 125 $ "Broyden_optim: Variable weight not implemented yet" 126 write(6,'(a,f8.4)') 127 $ "Broyden_optim: initial inverse jacobian: ", jinv0 128 write(6,*) 129 endif 130 131 call broyden_init(br,debug=broyden_debug) 132 initialization_done = .true. 133 134 endif 135 136 if ( frstme ) then 137 138 ! Find number of variables 139 ndeg = 0 140 do ia = 1,na 141 do n = 1,3 142 indi = 3*(ia-1) + n 143 if (VaryZmat(indi)) then 144 ndeg = ndeg + 1 145 endif 146 enddo 147 enddo 148 if ( varcel ) then 149 ndeg = ndeg + 6 ! Including stress 150 endif 151 if (Ionode) then 152 write(6,'(a,i6)') "Broyden_optim: No of elements: ", ndeg 153 endif 154 155 if ( varcel ) then 156 157 constant_volume = fdf_get("MD.ConstantVolume", .false.) 158 159 tarstr = fdf_block('MD.TargetStress',bfdf) 160 161 ! Look for target stress and read it if found, 162 ! otherwise generate it 163 if (tarstr) then 164 if (IOnode) then 165 write(6,'(/a,a)') 166 $ 'zm_broyden_optimizer: Reading %block MD.TargetStress', 167 . ' (units of MD.TargetPressure).' 168 endif 169 if (.not. fdf_bline(bfdf,pline)) 170 $ call die('zm_broyden_optimizer: ERROR in ' // 171 . 'MD.TargetStress block') 172 sxx = fdf_bvalues(pline,1) 173 syy = fdf_bvalues(pline,2) 174 szz = fdf_bvalues(pline,3) 175 sxy = fdf_bvalues(pline,4) 176 sxz = fdf_bvalues(pline,5) 177 syz = fdf_bvalues(pline,6) 178 call fdf_bclose(bfdf) 179 180 tstres(1,1) = - sxx * tp 181 tstres(2,2) = - syy * tp 182 tstres(3,3) = - szz * tp 183 tstres(1,2) = - sxy * tp 184 tstres(2,1) = - sxy * tp 185 tstres(1,3) = - sxz * tp 186 tstres(3,1) = - sxz * tp 187 tstres(2,3) = - syz * tp 188 tstres(3,2) = - syz * tp 189 else 190 write(6,'(/a,a)') 191 $ 'zm_broyden_optimizer: No target stress found, ', 192 . 'assuming hydrostatic MD.TargetPressure.' 193 do i = 1, 3 194 do j = 1, 3 195 tstres(i,j) = 0.0_dp 196 enddo 197 tstres(i,i) = - tp 198 enddo 199 endif 200 201 if (constant_volume) then 202 tstres(:,:) = 0.0_dp 203 if (IOnode) then 204 write(6,"(a)") "***Target stress set to zero " // 205 $ "for constant-volume calculation" 206 endif 207 endif 208 if (IOnode) then 209 write(6,"(/a)") 'cgvc_zmatrix: Target stress (kBar)' 210 do i=1, 3 211 write(6,"(a,2x,3f12.3)") 212 . 'cgvc_zmatrix:', (tstres(i,j)/kBar,j=1,3) 213 enddo 214 endif 215 216C Moduli of original cell vectors for fractional coor scaling back to au --- 217 do n = 1, 3 218 modcel(n) = 0.0_dp 219 do j = 1, 3 220 modcel(n) = modcel(n) + cell(j,n)*cell(j,n) 221 enddo 222 modcel(n) = dsqrt( modcel(n) ) 223 enddo 224 225C Scale factor for strain variables to share magnitude with coordinates 226C ---- (a length in Bohrs typical of bond lengths ..) 227 228 ! AG: This could better be volume^(1/3) by default 229 precon = fdf_get('MD.PreconditionVariableCell', 230 . 9.4486344d0,'Bohr') 231 232C Initialize absolute strain and save initial cell vectors 233C Initialization to 1 for numerical reasons, later substracted 234 235 strain(1:3,1:3) = 1.0_dp 236 cellin(1:3,1:3) = cell(1:3,1:3) 237 initial_volume = volcel(cellin) 238 239 endif ! varcel 240 241 frstme = .false. 242 endif ! First time 243 244C Allocate local memory 245 allocate(gfa(1:ndeg),gxa(1:ndeg)) 246 allocate(ftoln(1:ndeg),dxmaxn(1:ndeg)) 247 248! print *, "zm_broyden: cell in : ", cell 249! print *, "zm_broyden: strain in : ", strain 250 251 if ( varcel ) then 252 253 ! Inverse of matrix of cell vectors (transpose of) 254 call reclat( cell, celli, 0 ) 255 256C Transform coordinates and forces to fractional ---------------------------- 257C but scale them again to Bohr by using the (fix) moduli of the original ---- 258C lattice vectors (allows using maximum displacement as before) ------------- 259C convergence is checked here for input forces and compared with ftoln ------ 260 261 relaxd = .true. 262 ndegi = 0 263 ! Loop over degrees of freedom, scaling 264 ! only cartesian coordinates to fractional 265 do ia = 1,na 266 do n = 1,3 267 indi = 3*(ia-1) + n 268 if (VaryZmat(indi)) then 269 ndegi = ndegi + 1 270 if (ZmatType(indi).gt.1) then 271 ftoln(ndegi) = ZmatForceTolLen 272 dxmaxn(ndegi) = ZmatMaxDisplLen 273 else 274 ftoln(ndegi) = ZmatForceTolAng 275 dxmaxn(ndegi) = ZmatMaxDisplAng 276 endif 277 vi = iZmattoVars(indi) 278 if (vi.eq.0) then 279 force = ZmatForce(indi) 280 else 281 force = ZmatForceVar(vi) 282 endif 283 relaxd=relaxd.and.(dabs(force).lt.ftoln(ndegi)) 284 if (ZmatType(indi).gt.2) then 285C Cartesian coordinate 286 gxa(ndegi) = 0.0_dp 287 gfa(ndegi) = 0.0_dp 288 do i = 1,3 289 indi1 = 3*(ia-1)+i 290 gxa(ndegi) = gxa(ndegi)+ 291 . celli(i,n)*Zmat(indi1)*modcel(n) 292 if (i.eq.n) then 293 force1 = force 294 else 295 force1 = ZmatForce(indi1) 296 endif 297 gfa(ndegi) = gfa(ndegi)+ 298 . cell(i,n)*force1/modcel(n) 299 enddo 300 else 301 gxa(ndegi) = Zmat(indi) 302 gfa(ndegi) = force 303 endif 304 endif 305 enddo 306 enddo 307 308C Symmetrizing the stress tensor -------------------------------------------- 309 do i = 1,3 310 do j = i+1,3 311 stress(i,j) = 0.5_dp*( stress(i,j) + stress(j,i) ) 312 stress(j,i) = stress(i,j) 313 enddo 314 enddo 315 316C Subtract target stress 317 318 stress_dif = stress - tstres 319! 320! Take 1/3 of the trace out here if constant-volume needed 321! 322 if (constant_volume) then 323 trace = stress_dif(1,1) + stress_dif(2,2) + stress_dif(3,3) 324 do i=1,3 325 stress_dif(i,i) = stress_dif(i,i) - trace/3.0_dp 326 enddo 327 endif 328 329C Append stress (substracting target stress) and strain to gxa and gfa ------ 330C preconditioning: scale stress and strain so as to behave similarly to x,f - 331 do i = 1,3 332 do j = i,3 333 ndegi = ndegi + 1 334 gfa(ndegi) = -stress_dif(i,j)*volume/precon 335 gxa(ndegi) = strain(i,j) * precon 336 dxmaxn(ndegi) = ZmatMaxDisplLen 337 enddo 338 enddo 339 340C Check stress convergence 341 strtol = dabs(strtol) 342 do i = 1,3 343 do j = 1,3 344 relaxd = relaxd .and. 345 . ( dabs(stress_dif(i,j)) .lt. abs(strtol) ) 346 enddo 347 enddo 348 349 else ! FIXED CELL 350 351C Set number of degrees of freedom & variables 352 relaxd = .true. 353 ndegi = 0 354 do i = 1,na 355 do k = 1,3 356 indi = 3*(i-1)+k 357 if (VaryZmat(indi)) then 358 ndegi = ndegi + 1 359 gxa(ndegi) = Zmat(indi) 360 vi = iZmattoVars(indi) 361 if (vi.eq.0) then 362 force = ZmatForce(indi) 363 else 364 force = ZmatForceVar(vi) 365 endif 366 gfa(ndegi) = force 367 if (ZmatType(indi).eq.1) then 368 ftoln(ndegi) = ZmatForceTolAng 369 dxmaxn(ndegi) = ZmatMaxDisplAng 370 else 371 ftoln(ndegi) = ZmatForceTolLen 372 dxmaxn(ndegi) = ZmatMaxDisplLen 373 endif 374 relaxd=relaxd.and.(dabs(force).lt.ftoln(ndegi)) 375 endif 376 enddo 377 enddo 378 379 endif 380 381 if (relaxd) then 382 383 call local_dealloc() 384 385 return 386 end if 387 388 if (.not. broyden_is_setup(br)) then 389 call broyden_reset(br,ndeg,maxit,cycle_on_maxit, 390 $ jinv0,0.01_dp,dxmaxn) 391 endif 392 393 allocate(rnew(1:ndeg)) 394 395 weight = 1.0_dp 396 call broyden_step(br,gxa,gfa,w=weight,newx=rnew) 397 398 399 ! Transform back 400 if ( varcel ) then 401 402 ! New cell 403 indi = ndeg-6 404 do i = 1,3 405 do j = i,3 406 indi = indi + 1 407 strain(i,j) = rnew(indi) / precon 408 strain(j,i) = strain(i,j) 409 enddo 410 enddo 411 412 cell = cellin + matmul(strain-1.0_dp,cellin) 413 if (constant_volume) then 414 new_volume = volcel(cell) 415 if (Node.eq.0) write(6,"(a,f12.4)") 416 $ "Volume before coercion: ", new_volume/Ang**3 417 cell = cell * (initial_volume/new_volume)**(1.0_dp/3.0_dp) 418 endif 419 420C Output fractional coordinates to cartesian Bohr, and copy to xa ----------- 421 ndegi = 0 422 do ia=1,na 423 do k=1,3 424 indi = 3*(ia-1)+k 425 if (VaryZmat(indi)) then 426 ndegi = ndegi + 1 427 j = indi 428 do while (j.ne.0) 429 if (ZmatType(j).gt.2) then 430 Zmat(j) = 0.0_dp 431 do n=1,3 432 indi1 = 3*(ia-1)+n 433 ! Assume all three coords of this atom 434 ! are variables 435 ndi = ndegi + n - k 436 Zmat(j)=Zmat(j)+cell(k,n)*rnew(ndi)/modcel(n) 437 enddo 438 else 439 Zmat(j) = rnew(ndegi) 440 endif 441 Zmat(j) = Zmat(j)*coeffA(j) + coeffB(j) 442 j = iNextDept(j) 443 enddo 444 endif 445 enddo 446 enddo 447 call Zmat_to_Cartesian(xa) 448 449 else 450C Fixed cell 451C Return coordinates to correct arrays 452 ndegi = 0 453 do ia = 1,na 454 do k = 1,3 455 indi = 3*(ia-1)+k 456 if (VaryZmat(indi)) then 457 ndegi = ndegi + 1 458 j = indi 459 do while (j.ne.0) 460 Zmat(j) = rnew(ndegi)*coeffA(j)+coeffB(j) 461 j = iNextDept(j) 462 enddo 463 endif 464 enddo 465 enddo 466 call Zmat_to_Cartesian(xa) 467 endif 468 469! print *, "zm_broyden: cell out : ", cell 470! print *, "zm_broyden: strain out : ", strain 471 472 call local_dealloc() 473 474 contains 475 476 subroutine local_dealloc() 477 478C Deallocate local memory 479 if ( allocated(dxmaxn) ) deallocate(dxmaxn) 480 if ( allocated(ftoln) ) deallocate(ftoln) 481 if ( allocated(gxa) ) deallocate(gxa) 482 if ( allocated(gfa) ) deallocate(gfa) 483 if ( allocated(rnew) ) deallocate(rnew) 484 485 end subroutine local_dealloc 486 487 end subroutine zm_broyden_optimizer 488 489 end module m_zm_broyden_optim 490