8 |
|
#endif |
9 |
|
implicit none |
10 |
|
|
11 |
< |
contains |
11 |
> |
real(kind=dp), save :: rrf |
12 |
> |
real(kind=dp), save :: rt |
13 |
> |
real(kind=dp), save :: dielect |
14 |
> |
real(kind=dp), save :: rrfsq |
15 |
> |
real(kind=dp), save :: pre |
16 |
> |
logical, save :: rf_initialized = .false. |
17 |
|
|
18 |
+ |
contains |
19 |
+ |
|
20 |
+ |
subroutine initialize_rf() |
21 |
+ |
rrf = getRrf() |
22 |
+ |
rt = getRt() |
23 |
+ |
dielect = getDielect() |
24 |
+ |
|
25 |
+ |
rrfsq = rrf * rrf |
26 |
+ |
pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) |
27 |
+ |
|
28 |
+ |
rf_initialized = .true. |
29 |
+ |
|
30 |
+ |
return |
31 |
+ |
end subroutine initialize_rf |
32 |
+ |
|
33 |
|
subroutine accumulate_rf(atom1, atom2, rij, u_l) |
34 |
|
|
35 |
|
integer, intent(in) :: atom1, atom2 |
37 |
|
real (kind = dp), dimension(3, getNlocal()) :: u_l |
38 |
|
|
39 |
|
integer :: me1, me2 |
40 |
< |
real (kind = dp) :: rrf, rt, taper, mu1, mu2 |
40 |
> |
real (kind = dp) :: taper, mu1, mu2 |
41 |
|
real (kind = dp), dimension(3) :: ul1 |
42 |
< |
real (kind = dp), dimension(3) :: ul2 |
42 |
> |
real (kind = dp), dimension(3) :: ul2 |
43 |
> |
|
44 |
> |
if (.not.rf_initialized) then |
45 |
> |
write(default_error,*) 'Reaction field not initialized!' |
46 |
> |
return |
47 |
> |
endif |
48 |
|
|
24 |
– |
rrf = getRrf() |
25 |
– |
|
49 |
|
if (rij.le.rrf) then |
50 |
< |
|
28 |
< |
rt = getRt() |
29 |
< |
|
50 |
> |
|
51 |
|
if (rij.lt.rt) then |
52 |
|
taper = 1.0d0 |
53 |
|
else |
116 |
|
return |
117 |
|
end subroutine accumulate_self_rf |
118 |
|
|
119 |
< |
subroutine reaction_field(a1, mu1, u_l, rfpot, t, do_pot) |
119 |
> |
subroutine reaction_field_final(a1, mu1, u_l, rfpot, t, do_pot) |
120 |
|
|
121 |
+ |
integer, intent(in) :: a1 |
122 |
+ |
real (kind=dp), intent(in) :: mu1 |
123 |
+ |
real (kind=dp), intent(inout) :: rfpot |
124 |
+ |
logical, intent(in) :: do_pot |
125 |
+ |
real (kind = dp), dimension(3, getNlocal()) :: u_l |
126 |
+ |
real (kind = dp), dimension(3, getNlocal()) :: t |
127 |
+ |
|
128 |
+ |
if (.not.rf_initialized) then |
129 |
+ |
write(default_error,*) 'Reaction field not initialized!' |
130 |
+ |
return |
131 |
+ |
endif |
132 |
+ |
|
133 |
|
! compute torques on dipoles: |
134 |
|
! pre converts from mu in units of debye to kcal/mol |
135 |
|
|
103 |
– |
rrf = getRrf() |
104 |
– |
dielect = getDielect() |
105 |
– |
rrfsq = rrf * rrf |
106 |
– |
pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) |
107 |
– |
|
136 |
|
! The torque contribution is dipole cross reaction_field |
137 |
|
|
138 |
|
t(1,a1) = t(1,a1) + pre*mu1*(u_l(2,a1)*rf(3,a1) - u_l(3,a1)*rf(2,a1)) |
147 |
|
endif |
148 |
|
|
149 |
|
return |
150 |
< |
end subroutine reaction_field |
150 |
> |
end subroutine reaction_field_final |
151 |
|
|
152 |
|
subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, f, do_stress) |
153 |
|
|
160 |
|
|
161 |
|
real (kind = dp), dimension(3) :: ul1 |
162 |
|
real (kind = dp), dimension(3) :: ul2 |
163 |
< |
real (kind = dp) :: dtdr, rrfsq, prerf |
163 |
> |
real (kind = dp) :: dtdr |
164 |
|
real (kind = dp) :: dudx, dudy, dudz, u1dotu2 |
165 |
+ |
integer :: me1, me2 |
166 |
+ |
real (kind = dp) :: mu1, mu2 |
167 |
|
|
168 |
< |
rrf = getRrf() |
169 |
< |
|
168 |
> |
if (.not.rf_initialized) then |
169 |
> |
write(default_error,*) 'Reaction field not initialized!' |
170 |
> |
return |
171 |
> |
endif |
172 |
> |
|
173 |
|
if (rij.le.rrf) then |
174 |
|
|
142 |
– |
rrfsq = rrf * rrf |
143 |
– |
dielect = getDielect() |
144 |
– |
prerf = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) |
145 |
– |
rt = getRt() |
146 |
– |
|
175 |
|
if (rij.lt.rt) then |
176 |
|
dtdr = 0.0d0 |
177 |
|
else |
205 |
|
|
206 |
|
u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3) |
207 |
|
|
208 |
< |
dudx = - prerf*mu1*mu2*u1dotu2*dtdr*d(1)/rij |
209 |
< |
dudy = - prerf*mu1*mu2*u1dotu2*dtdr*d(2)/rij |
210 |
< |
dudz = - prerf*mu1*mu2*u1dotu2*dtdr*d(3)/rij |
208 |
> |
dudx = - pre*mu1*mu2*u1dotu2*dtdr*d(1)/rij |
209 |
> |
dudy = - pre*mu1*mu2*u1dotu2*dtdr*d(2)/rij |
210 |
> |
dudz = - pre*mu1*mu2*u1dotu2*dtdr*d(3)/rij |
211 |
|
|
212 |
|
#ifdef IS_MPI |
213 |
|
f_Row(1,atom1) = f_Row(1,atom1) + dudx |
238 |
|
tau_Temp(8) = tau_Temp(8) + dudz * d(2) |
239 |
|
tau_Temp(9) = tau_Temp(9) + dudz * d(3) |
240 |
|
virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) |
241 |
< |
endif |
214 |
< |
|
241 |
> |
endif |
242 |
|
endif |
243 |
|
|
244 |
|
return |