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