1! This file is part of xtb. 2! 3! Copyright (C) 2019-2020 Stefan Grimme 4! 5! xtb is free software: you can redistribute it and/or modify it under 6! the terms of the GNU Lesser General Public License as published by 7! the Free Software Foundation, either version 3 of the License, or 8! (at your option) any later version. 9! 10! xtb is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU Lesser General Public License for more details. 14! 15! You should have received a copy of the GNU Lesser General Public License 16! along with xtb. If not, see <https://www.gnu.org/licenses/>. 17 18 module xtb_gfnff_shake 19 implicit none 20 integer :: shake_mode 21 integer :: ncons = 0 22 integer :: nconsu = 0 23 integer, allocatable :: conslist(:,:) 24 real*8, allocatable :: distcons(:) 25 real*8, allocatable :: dro(:,:) 26 real*8, allocatable :: dr (:,:) 27 integer, parameter :: maxcyc = 5000 28 real*8, parameter :: tolshake = 1.d-6 29 30 contains 31 32 subroutine init_shake(nat,at,xyz,topo) 33 use xtb_gfnff_topology, only : TGFFTopology 34 implicit none 35 type(TGFFTopology), intent(in) :: topo 36 integer :: nat,at(nat) 37 real*8 :: xyz(3,nat) 38 real*8 :: wbo(nat,nat) 39 integer :: iat,jat,i,j,jmin 40 real*8 :: minrij,rij,wthr 41 real*8 :: drij(3) 42 43 integer list(nat*(nat+1)/2),lin,ij 44 45 ncons = topo%nbond 46 allocate(conslist(2,ncons),distcons(ncons),& 47 & dro(3,ncons),dr(4,ncons)) 48 49 conslist(1:2,1:ncons)=topo%blist(1:2,1:ncons) 50 do i = 1, ncons 51 iat = conslist(1,i) 52 jat = conslist(2,i) 53 drij(1:3) = xyz(1:3,iat) - xyz(1:3,jat) 54 distcons(i) = drij(1)**2 + drij(2)**2 + drij(3)**2 55 enddo 56 57 return 58 end subroutine init_shake 59 60 subroutine do_shake(nat,xyzo,xyz,velo,acc,mass,tstep) 61 implicit none 62 integer nat 63 real*8 xyzo(3,nat),xyz(3,nat),velo(3,nat),acc(3,nat),mass(nat) 64 real*8 virsh(3),tstep 65 real*8 xyzt(3,nat) 66 real*8 vel(3) 67 integer i,icyc 68 integer iat,jat 69 logical conv 70 real*8 maxdev 71 real*8 r, dev,dist,denom 72 real*8 gcons, rmi, rmj 73 real*8 tau1,tau2 74 integer jmaxdev 75 76 conv = .false. 77 icyc = 0 78 79 xyzt = xyz 80 81 tau1 = 1.d0/tstep 82 tau2 = tau1*tau1 83 84 do i = 1, ncons 85 iat = conslist(1,i) 86 jat = conslist(2,i) 87 dro(1:3,i) = xyzo(1:3,iat) - xyzo(1:3,jat) 88 enddo 89 90!100 continue 91 do while (.not.conv.and.icyc.le.maxcyc) !goto 100 92 93 maxdev = 0.d0 94 95 do i = 1, ncons 96 iat = conslist(1,i) 97 jat = conslist(2,i) 98 dr(1:3,i) = xyzt(1:3,iat) - xyzt(1:3,jat) 99 dr(4,i) = dr(1,i)**2 + dr(2,i)**2 + dr(3,i)**2 100 dist = distcons(i) 101 dev = abs(dr(4,i) - dist) / dist 102 if(dev.gt.maxdev) then 103 maxdev = dev 104 jmaxdev = i 105 endif 106 enddo 107 108 if(maxdev.lt.tolshake) conv = .true. 109 110 if(.not.conv) then 111 112 do i = 1, ncons 113 114 iat = conslist(1,i) 115 jat = conslist(2,i) 116 dist = distcons(i) 117 118 rmi = 1.d0/mass(iat) 119 rmj = 1.d0/mass(jat) 120 121 denom = 2.d0*(rmi+rmj)*(dr(1,i)*dro(1,i)+dr(2,i)*dro(2,i)+& 122 & dr(3,i)*dro(3,i)) 123 124 gcons = (dist-dr(4,i))/denom 125 126 xyzt(1:3,iat) = xyzt(1:3,iat) + rmi*gcons*dro(1:3,i) 127 xyzt(1:3,jat) = xyzt(1:3,jat) - rmj*gcons*dro(1:3,i) 128 129 enddo 130 131 endif 132 133 icyc = icyc + 1 134 135 end do 136 !if(.not.conv.and.icyc.le.maxcyc) goto 100 137 138 if(conv) then 139 velo = velo + (xyzt-xyz)*tau1 140 acc = acc + (xyzt-xyz)*tau2 141 xyz = xyzt 142 else 143 write(*,*)'SHAKE did not converge! maxdev=',maxdev 144 !if(maxdev.gt.1.d-3) stop 'SHAKE error too large' 145 endif 146 147 return 148 end subroutine do_shake 149 150 end module xtb_gfnff_shake 151