-
Notifications
You must be signed in to change notification settings - Fork 1
/
egrad_dg_evb.f90
210 lines (196 loc) · 6.6 KB
/
egrad_dg_evb.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! CARACAL - Ring polymer molecular dynamics and rate constant calculations
! on black-box generated potential energy surfaces
!
! Copyright (c) 2023 by Julien Steffen ([email protected])
! Stefan Grimme ([email protected]) (QMDFF code)
!
! Permission is hereby granted, free of charge, to any person obtaining a
! copy of this software and associated documentation files (the "Software"),
! to deal in the Software without restriction, including without limitation
! the rights to use, copy, modify, merge, publish, distribute, sublicense,
! and/or sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! subroutine egrad_dg_evb: calculate energy and gradient of the Distributed
! Gaussian (DG)-EVB coupling term
!
! part of EVB
!
subroutine egrad_dg_evb(xyz2,e_qmdff1,e_qmdff2,ediff,e_evb,g_evb)
use general
use evb_mod
implicit none
! the actual structure
real(kind=8), intent(inout) :: xyz2(3,natoms)
! energies of both QMDFFs and their difference
real(kind=8), intent(inout) :: e_qmdff1,e_qmdff2,ediff
! the calculated DG-EVB energy and gradient
real(kind=8), intent(out) :: e_evb,g_evb(3,natoms)
! loop indices
integer::i,j,k
! the DG-EVB mode
integer::dg_evb_mode
! size of the DG-EVB coefficient matrix
integer::mat_size
! the current internal coordinates
real(kind=8),dimension(nat6)::int_path
! general EVB calculation parameters
real(kind=8)::offdiag,off4,root2,deldiscr,delsqrt,root
! the EVB coupling strength and its internal and cartesian derivatives
real(kind=8)::V12,g_V12_int(nat6),g_V12(3,natoms)
! energies of the QMDFFs without shifting
real(kind=8)::e,e_two
! for numrical gradient
real(kind=8)::e_lower,e_upper,step
! the internal gradient
real(kind=8)::g_int(nat6)
! avoid imaginary coupling values
logical::unset_V12
! calculate the gradient (analytical and mumerical):
!
dg_evb_mode=dg_mode
if (dg_evb_mode .eq. 1) then
mat_size=dg_evb_points
else if (dg_evb_mode .eq. 2) then
mat_size=dg_evb_points*(1+nat6)
else if (dg_evb_mode .eq. 3) then
mat_size=dg_evb_points*(1+nat6+(nat6)*(nat6+1)/2)
end if
if (num_grad .eqv. .true.) then
!
! The numerical DG-EVB gradient
!
step=num_grad_step
do i=1,natoms
do j=1,3
do k=1,2
if (k.eq.1) then
xyz2(j,i)=xyz2(j,i)+step
else if (k.eq.2) then
xyz2(j,i)=xyz2(j,i)-2d0*step
end if
call ff_eg(n_one,at,xyz2,e,g_one)
call ff_nonb(n_one,at,xyz2,q,r0ab,zab,r094_mod,&
&sr42,c6xy,e,g_one)
call ff_hb(n_one,at,xyz2,e,g_one)
call ff_eg_two(n_one,at,xyz2,e_two,g_two)
call ff_nonb_two(n_one,at,xyz2,q_two,r0ab,zab,&
&r094_mod,sr42,c6xy_two,e_two,g_two)
call ff_hb_two(n_one,at,xyz2,e_two,g_two)
call xyz_2int(xyz2,int_path,natoms)
e_qmdff1=e+E_zero1
e_qmdff2=e_two+E_zero2
V12=0
!
! call subroutine for calculating the DG-EVB-resonance integral
!
call sum_v12(dg_evb_mode,mat_size,alph_opt,int_path,V12)
!
! prohibit negative coupling-element values
!
root=(0.5d0*(e_qmdff1-e_qmdff2))*(0.5d0*(e_qmdff1-e_qmdff2))+V12
if (root .le. 0) then
e=0.5d0*(e_qmdff1+e_qmdff2)
else
e=0.5d0*(e_qmdff1+e_qmdff2)-sqrt(root)
end if
if (k.eq.1) then
e_upper=e
else if (k.eq.2) then
e_lower=e
end if
end do
xyz2(j,i)=xyz2(j,i)+step
g_evb(j,i)=(e_upper-e_lower)/(2d0*step)
g_evb(j,i)=g_evb(j,i)
end do
end do
else
!
! The analytical DG-EVB gradient
! Newly implemented: 22.07.2018!
! Now the calculation is done complete analytically, like in
! the case of RP-EVB
call xyz_2int(xyz2,int_path,natoms)
!
! calculate the coupling strengh
!
call sum_v12(dg_evb_mode,mat_size,alph_opt,int_path,V12)
!
! calculate the coupling gradient in internal coordinates!
!
call sum_dv12(dg_evb_mode,mat_size,alph_opt,int_path,g_V12_int)
!
! Convert the squared coupling gradient to cartesian coordinates
!
call int2grad(xyz2,int_path,g_V12_int,g_V12)
!
! Now, calculate the EVB-QMDFF gradient using the well known formula
! Using some mathematical tricks, we can avoid the root of the coupling term!
!
offdiag=V12
off4=4.d0*offdiag
!
! If the root argument becomes negative, set the coupling to zero..
!
unset_V12=.false.
if (ediff*ediff+off4 .lt. 0.d0) unset_V12=.true.
if (unset_V12 .eqv..false.) then
root2=sqrt(ediff*ediff+off4)
end if
do i=1,n_one
do j=1,3
deldiscr=ediff*(g_one(j,i)-g_two(j,i))+2.d0*g_V12(j,i)
if (unset_V12 .eqv. .false.) then
delsqrt=deldiscr/root2
else
delsqrt=0.d0
end if
g_evb(j,i)=0.5d0*(g_one(j,i)+g_two(j,i)-delsqrt)
end do
end do
end if
!
! Optional test: write components of internal gradients to file
!
if (int_grad_plot) then
g_int=0.d0
call grad2int(xyz2,int_path,g_int,g_evb)
write(192,*) sum(g_int**2)
end if
!
! calculate energy of the undistorted structure
!
call xyz_2int(xyz2,int_path,natoms)
V12=0
!
! Write internal coordinates to file if option is activated
! Print a third column for splot an convert them into Angstrom
!
if (int_coord_plot) then
write(191,*) int_path*bohr,0.0d0
end if
call sum_v12(dg_evb_mode,mat_size,alph_opt,int_path,V12)
root=(0.5d0*(e_qmdff1-e_qmdff2))*(0.5d0*(e_qmdff1-e_qmdff2))+V12
if (root .le. 0) then
e_evb=0.5d0*(e_qmdff1+e_qmdff2)
else
e_evb=0.5d0*(e_qmdff1+e_qmdff2)-sqrt(root)
end if
return
end subroutine egrad_dg_evb