46 |
|
use vector_class |
47 |
|
use simulation |
48 |
|
use status |
49 |
+ |
use electrostatic_module |
50 |
|
#ifdef IS_MPI |
51 |
|
use mpiSimulation |
52 |
|
#endif |
59 |
|
real(kind=dp), save :: dielect = 1.0_dp |
60 |
|
real(kind=dp), save :: rrfsq = 1.0_dp |
61 |
|
real(kind=dp), save :: pre |
62 |
+ |
|
63 |
|
logical, save :: haveCutoffs = .false. |
62 |
– |
logical, save :: haveMomentMap = .false. |
64 |
|
logical, save :: haveDielectric = .false. |
65 |
|
|
65 |
– |
type :: MomentList |
66 |
– |
real(kind=DP) :: dipole_moment = 0.0_DP |
67 |
– |
end type MomentList |
68 |
– |
|
69 |
– |
type(MomentList), dimension(:),allocatable :: MomentMap |
70 |
– |
|
66 |
|
PUBLIC::initialize_rf |
67 |
|
PUBLIC::setCutoffsRF |
68 |
|
PUBLIC::accumulate_rf |
77 |
|
|
78 |
|
dielect = this_dielect |
79 |
|
|
80 |
< |
pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) |
80 |
> |
pre = pre22 * 2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) |
81 |
|
|
82 |
|
haveDielectric = .true. |
83 |
|
|
92 |
|
rt = this_rt |
93 |
|
|
94 |
|
rrfsq = rrf * rrf |
95 |
< |
pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) |
95 |
> |
pre = pre22 * 2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) |
96 |
|
|
97 |
|
haveCutoffs = .true. |
98 |
|
|
99 |
|
end subroutine setCutoffsRF |
105 |
– |
|
106 |
– |
subroutine createMomentMap(status) |
107 |
– |
integer :: nAtypes |
108 |
– |
integer :: status |
109 |
– |
integer :: i |
110 |
– |
real (kind=DP) :: thisDP |
111 |
– |
logical :: thisProperty |
112 |
– |
|
113 |
– |
status = 0 |
114 |
– |
|
115 |
– |
nAtypes = getSize(atypes) |
116 |
– |
|
117 |
– |
if (nAtypes == 0) then |
118 |
– |
status = -1 |
119 |
– |
return |
120 |
– |
end if |
121 |
– |
|
122 |
– |
if (.not. allocated(MomentMap)) then |
123 |
– |
allocate(MomentMap(nAtypes)) |
124 |
– |
endif |
125 |
– |
|
126 |
– |
!!do i = 1, nAtypes |
127 |
– |
!! |
128 |
– |
!! call getElementProperty(atypes, i, "is_Dipole", thisProperty) |
129 |
– |
!! |
130 |
– |
!! if (thisProperty) then |
131 |
– |
!! call getElementProperty(atypes, i, "dipole_moment", thisDP) |
132 |
– |
!! MomentMap(i)%dipole_moment = thisDP |
133 |
– |
!! endif |
134 |
– |
!! |
135 |
– |
!!end do |
136 |
– |
|
137 |
– |
haveMomentMap = .true. |
138 |
– |
|
139 |
– |
end subroutine createMomentMap |
100 |
|
|
101 |
|
subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper) |
102 |
|
|
115 |
|
if ((.not.haveDielectric).or.(.not.haveCutoffs)) then |
116 |
|
write(default_error,*) 'Reaction field not initialized!' |
117 |
|
return |
118 |
< |
endif |
159 |
< |
|
160 |
< |
if (.not.haveMomentMap) then |
161 |
< |
localError = 0 |
162 |
< |
call createMomentMap(localError) |
163 |
< |
if ( localError .ne. 0 ) then |
164 |
< |
call handleError("reaction-field", "MomentMap creation failed!") |
165 |
< |
return |
166 |
< |
end if |
167 |
< |
endif |
168 |
< |
|
118 |
> |
endif |
119 |
|
|
120 |
|
#ifdef IS_MPI |
121 |
|
me1 = atid_Row(atom1) |
139 |
|
ul2(3) = eFrame(9,atom2) |
140 |
|
#endif |
141 |
|
|
142 |
< |
mu1 = MomentMap(me1)%dipole_moment |
143 |
< |
mu2 = MomentMap(me2)%dipole_moment |
142 |
> |
mu1 = getDipoleMoment(me1) |
143 |
> |
mu2 = getDipoleMoment(me2) |
144 |
|
|
145 |
|
#ifdef IS_MPI |
146 |
|
rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper |
194 |
|
return |
195 |
|
endif |
196 |
|
|
247 |
– |
if (.not.haveMomentMap) then |
248 |
– |
localError = 0 |
249 |
– |
call createMomentMap(localError) |
250 |
– |
if ( localError .ne. 0 ) then |
251 |
– |
call handleError("reaction-field", "MomentMap creation failed!") |
252 |
– |
return |
253 |
– |
end if |
254 |
– |
endif |
255 |
– |
|
197 |
|
! compute torques on dipoles: |
198 |
|
! pre converts from mu in units of debye to kcal/mol |
199 |
|
|
236 |
|
return |
237 |
|
endif |
238 |
|
|
298 |
– |
if (.not.haveMomentMap) then |
299 |
– |
localError = 0 |
300 |
– |
call createMomentMap(localError) |
301 |
– |
if ( localError .ne. 0 ) then |
302 |
– |
call handleError("reaction-field", "MomentMap creation failed!") |
303 |
– |
return |
304 |
– |
end if |
305 |
– |
endif |
306 |
– |
|
239 |
|
if (rij.le.rrf) then |
240 |
|
|
241 |
|
if (rij.lt.rt) then |
267 |
|
ul2(3) = eFrame(9,atom2) |
268 |
|
#endif |
269 |
|
|
270 |
< |
mu1 = MomentMap(me1)%dipole_moment |
271 |
< |
mu2 = MomentMap(me2)%dipole_moment |
270 |
> |
mu1 = getDipoleMoment(me1) |
271 |
> |
mu2 = getDipoleMoment(me2) |
272 |
|
|
273 |
|
u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3) |
274 |
|
|