ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/electrostatic.F90 (file contents):
Revision 2339 by chrisfen, Wed Sep 28 18:47:17 2005 UTC vs.
Revision 2390 by chrisfen, Wed Oct 19 19:24:40 2005 UTC

# Line 54 | Line 54 | module electrostatic_module
54  
55    PRIVATE
56  
57 +
58   #define __FORTRAN90
59 + #include "UseTheForce/DarkSide/fInteractionMap.h"
60   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61  
62 +
63    !! these prefactors convert the multipole interactions into kcal / mol
64    !! all were computed assuming distances are measured in angstroms
65    !! Charge-Charge, assuming charges are measured in electrons
# Line 75 | Line 78 | module electrostatic_module
78    integer, save :: summationMethod = NONE
79    logical, save :: summationMethodChecked = .false.
80    real(kind=DP), save :: defaultCutoff = 0.0_DP
81 +  real(kind=DP), save :: defaultCutoff2 = 0.0_DP
82    logical, save :: haveDefaultCutoff = .false.
83    real(kind=DP), save :: dampingAlpha = 0.0_DP
84    logical, save :: haveDampingAlpha = .false.
85 <  real(kind=DP), save :: dielectric = 0.0_DP
85 >  real(kind=DP), save :: dielectric = 1.0_DP
86    logical, save :: haveDielectric = .false.
87    real(kind=DP), save :: constERFC = 0.0_DP
88    real(kind=DP), save :: constEXP = 0.0_DP
89    logical, save :: haveDWAconstants = .false.
90 <  real(kind=dp), save :: rcuti = 0.0_dp
91 <  real(kind=dp), save :: rcuti2 = 0.0_dp
92 <  real(kind=dp), save :: rcuti3 = 0.0_dp
93 <  real(kind=dp), save :: rcuti4 = 0.0_dp
94 <  real(kind=dp), save :: alphaPi = 0.0_dp
95 <  real(kind=dp), save :: invRootPi = 0.0_dp
96 <  
90 >  real(kind=dp), save :: rcuti = 0.0_DP
91 >  real(kind=dp), save :: rcuti2 = 0.0_DP
92 >  real(kind=dp), save :: rcuti3 = 0.0_DP
93 >  real(kind=dp), save :: rcuti4 = 0.0_DP
94 >  real(kind=dp), save :: alphaPi = 0.0_DP
95 >  real(kind=dp), save :: invRootPi = 0.0_DP
96 >  real(kind=dp), save :: rrf = 1.0_DP
97 >  real(kind=dp), save :: rt = 1.0_DP
98 >  real(kind=dp), save :: rrfsq = 1.0_DP
99 >  real(kind=dp), save :: preRF = 0.0_DP
100 >  logical, save :: preRFCalculated = .false.
101 >
102   #ifdef __IFC
103   ! error function for ifc version > 7.
104    double precision, external :: derfc
# Line 99 | Line 108 | module electrostatic_module
108    public :: setElectrostaticCutoffRadius
109    public :: setDampedWolfAlpha
110    public :: setReactionFieldDielectric
111 +  public :: setReactionFieldPrefactor
112    public :: newElectrostaticType
113    public :: setCharge
114    public :: setDipoleMoment
# Line 109 | Line 119 | module electrostatic_module
119    public :: getDipoleMoment
120    public :: pre22
121    public :: destroyElectrostaticTypes
122 +  public :: accumulate_rf
123 +  public :: accumulate_self_rf
124 +  public :: reaction_field_final
125 +  public :: rf_correct_forces
126  
127    type :: Electrostatic
128       integer :: c_ident
# Line 138 | Line 152 | contains
152  
153    end subroutine setElectrostaticSummationMethod
154  
155 <  subroutine setElectrostaticCutoffRadius(thisRcut)
155 >  subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
156      real(kind=dp), intent(in) :: thisRcut
157 +    real(kind=dp), intent(in) :: thisRsw
158      defaultCutoff = thisRcut
159 +    rrf = defaultCutoff
160 +    rt = thisRsw
161      haveDefaultCutoff = .true.
162    end subroutine setElectrostaticCutoffRadius
163  
# Line 155 | Line 172 | contains
172      dielectric = thisDielectric
173      haveDielectric = .true.
174    end subroutine setReactionFieldDielectric
175 +
176 +  subroutine setReactionFieldPrefactor
177 +    if (haveDefaultCutoff .and. haveDielectric) then
178 +       defaultCutoff2 = defaultCutoff*defaultCutoff
179 +       preRF = pre22 * 2.0d0*(dielectric-1.0d0) / &
180 +            ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
181 +       preRFCalculated = .true.
182 +    else if (.not.haveDefaultCutoff) then
183 +       call handleError("setReactionFieldPrefactor", "Default cutoff not set")
184 +    else
185 +       call handleError("setReactionFieldPrefactor", "Dielectric not set")
186 +    endif
187 +  end subroutine setReactionFieldPrefactor
188  
189    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
190         is_SplitDipole, is_Quadrupole, is_Tap, status)
# Line 392 | Line 422 | contains
422            constERFC = derfc(dampingAlpha*defaultCutoff)
423            invRootPi = 0.56418958354775628695d0
424            alphaPi = 2*dampingAlpha*invRootPi
425 <          
425 >  
426            haveDWAconstants = .true.
427         endif
428      endif
# Line 448 | Line 478 | contains
478      real (kind=dp) :: dudx, dudy, dudz
479      real (kind=dp) :: scale, sc2, bigR
480      real (kind=dp) :: varERFC, varEXP
481 +    real (kind=dp) :: limScale
482  
483      if (.not.allocated(ElectrostaticMap)) then
484         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 471 | Line 502 | contains
502      !! some variables we'll need independent of electrostatic type:
503  
504      riji = 1.0d0 / rij
505 <
505 >  
506      xhat = d(1) * riji
507      yhat = d(2) * riji
508      zhat = d(3) * riji
# Line 617 | Line 648 | contains
648               vpair = vpair + vterm
649               epot = epot + sw*vterm
650              
651 <             dudr  = -sw*pre11*q_i*q_j * (riji*riji*riji - rcuti2*rcuti)
651 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
652              
653               dudx = dudx + dudr * d(1)
654               dudy = dudy + dudr * d(2)
# Line 631 | Line 662 | contains
662               vpair = vpair + vterm
663               epot = epot + sw*vterm
664              
665 <             dudr  = -sw*pre11*q_i*q_j * ( riji*(varERFC*riji*riji &
666 <                                                 + alphaPi*varEXP) &
667 <                                         - rcuti*(constERFC*rcuti2 &
668 <                                                 + alphaPi*constEXP) )
665 >             dudr  = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
666 >                                                  + alphaPi*varEXP) &
667 >                                                 - (constERFC*rcuti2 &
668 >                                                    + alphaPi*constEXP)) )
669              
670               dudx = dudx + dudr * d(1)
671               dudy = dudy + dudr * d(2)
# Line 1080 | Line 1111 | contains
1111  
1112      if (do_pot) then
1113   #ifdef IS_MPI
1114 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1115 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1114 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1115 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1116   #else
1117         pot = pot + epot
1118   #endif
# Line 1186 | Line 1217 | contains
1217      return
1218    end subroutine doElectrostaticPair
1219  
1189 –  !! calculates the switching functions and their derivatives for a given
1190 –  subroutine calc_switch(r, mu, scale, dscale)
1191 –
1192 –    real (kind=dp), intent(in) :: r, mu
1193 –    real (kind=dp), intent(inout) :: scale, dscale
1194 –    real (kind=dp) :: rl, ru, mulow, minRatio, temp, scaleVal
1195 –
1196 –    ! distances must be in angstroms
1197 –    rl = 2.75d0
1198 –    ru = 3.75d0
1199 –    mulow = 0.0d0 !3.3856d0 ! 1.84 * 1.84
1200 –    minRatio = mulow / (mu*mu)
1201 –    scaleVal = 1.0d0 - minRatio
1202 –    
1203 –    if (r.lt.rl) then
1204 –       scale = minRatio
1205 –       dscale = 0.0d0
1206 –    elseif (r.gt.ru) then
1207 –       scale = 1.0d0
1208 –       dscale = 0.0d0
1209 –    else
1210 –       scale = 1.0d0 - scaleVal*((ru + 2.0d0*r - 3.0d0*rl) * (ru-r)**2) &
1211 –                        / ((ru - rl)**3)
1212 –       dscale = -scaleVal * 6.0d0 * (r-ru)*(r-rl)/((ru - rl)**3)    
1213 –    endif
1214 –        
1215 –    return
1216 –  end subroutine calc_switch
1217 –
1220    subroutine destroyElectrostaticTypes()
1221  
1222      if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1223  
1224    end subroutine destroyElectrostaticTypes
1225  
1226 +  subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper)
1227 +
1228 +    integer, intent(in) :: atom1, atom2
1229 +    real (kind = dp), intent(in) :: rij
1230 +    real (kind = dp), dimension(9,nLocal) :: eFrame
1231 +
1232 +    integer :: me1, me2
1233 +    real (kind = dp), intent(in) :: taper
1234 +    real (kind = dp):: mu1, mu2
1235 +    real (kind = dp), dimension(3) :: ul1
1236 +    real (kind = dp), dimension(3) :: ul2  
1237 +
1238 +    integer :: localError
1239 +
1240 + #ifdef IS_MPI
1241 +    me1 = atid_Row(atom1)
1242 +    ul1(1) = eFrame_Row(3,atom1)
1243 +    ul1(2) = eFrame_Row(6,atom1)
1244 +    ul1(3) = eFrame_Row(9,atom1)
1245 +
1246 +    me2 = atid_Col(atom2)
1247 +    ul2(1) = eFrame_Col(3,atom2)
1248 +    ul2(2) = eFrame_Col(6,atom2)
1249 +    ul2(3) = eFrame_Col(9,atom2)
1250 + #else
1251 +    me1 = atid(atom1)
1252 +    ul1(1) = eFrame(3,atom1)
1253 +    ul1(2) = eFrame(6,atom1)
1254 +    ul1(3) = eFrame(9,atom1)
1255 +
1256 +    me2 = atid(atom2)
1257 +    ul2(1) = eFrame(3,atom2)
1258 +    ul2(2) = eFrame(6,atom2)
1259 +    ul2(3) = eFrame(9,atom2)
1260 + #endif
1261 +
1262 +    mu1 = getDipoleMoment(me1)
1263 +    mu2 = getDipoleMoment(me2)
1264 +
1265 + #ifdef IS_MPI
1266 +    rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper
1267 +    rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper
1268 +    rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper
1269 +
1270 +    rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu1*taper
1271 +    rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper
1272 +    rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*taper
1273 + #else
1274 +    rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper
1275 +    rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper
1276 +    rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper
1277 +
1278 +    rf(1,atom2) = rf(1,atom2) + ul1(1)*mu1*taper
1279 +    rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper
1280 +    rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper    
1281 + #endif
1282 +    return  
1283 +  end subroutine accumulate_rf
1284 +
1285 +  subroutine accumulate_self_rf(atom1, mu1, eFrame)
1286 +
1287 +    integer, intent(in) :: atom1
1288 +    real(kind=dp), intent(in) :: mu1
1289 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1290 +
1291 +    !! should work for both MPI and non-MPI version since this is not pairwise.
1292 +    rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
1293 +    rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
1294 +    rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1
1295 +
1296 +    return
1297 +  end subroutine accumulate_self_rf
1298 +
1299 +  subroutine reaction_field_final(a1, mu1, eFrame, rfpot, t, do_pot)
1300 +
1301 +    integer, intent(in) :: a1
1302 +    real (kind=dp), intent(in) :: mu1
1303 +    real (kind=dp), intent(inout) :: rfpot
1304 +    logical, intent(in) :: do_pot
1305 +    real (kind = dp), dimension(9,nLocal) :: eFrame
1306 +    real (kind = dp), dimension(3,nLocal) :: t
1307 +
1308 +    integer :: localError
1309 +
1310 +    if (.not.preRFCalculated) then
1311 +       call setReactionFieldPrefactor()
1312 +    endif
1313 +
1314 +    ! compute torques on dipoles:
1315 +    ! pre converts from mu in units of debye to kcal/mol
1316 +
1317 +    ! The torque contribution is dipole cross reaction_field  
1318 +
1319 +    t(1,a1) = t(1,a1) + preRF*mu1*(eFrame(6,a1)*rf(3,a1) - &
1320 +                                   eFrame(9,a1)*rf(2,a1))
1321 +    t(2,a1) = t(2,a1) + preRF*mu1*(eFrame(9,a1)*rf(1,a1) - &
1322 +                                   eFrame(3,a1)*rf(3,a1))
1323 +    t(3,a1) = t(3,a1) + preRF*mu1*(eFrame(3,a1)*rf(2,a1) - &
1324 +                                   eFrame(6,a1)*rf(1,a1))
1325 +
1326 +    ! the potential contribution is -1/2 dipole dot reaction_field
1327 +
1328 +    if (do_pot) then
1329 +       rfpot = rfpot - 0.5d0 * preRF * mu1 * &
1330 +            (rf(1,a1)*eFrame(3,a1) + rf(2,a1)*eFrame(6,a1) + &
1331 +             rf(3,a1)*eFrame(9,a1))
1332 +    endif
1333 +
1334 +    return
1335 +  end subroutine reaction_field_final
1336 +
1337 +  subroutine rf_correct_forces(atom1, atom2, d, rij, eFrame, taper, f, fpair)
1338 +
1339 +    integer, intent(in) :: atom1, atom2
1340 +    real(kind=dp), dimension(3), intent(in) :: d
1341 +    real(kind=dp), intent(in) :: rij, taper
1342 +    real( kind = dp ), dimension(9,nLocal) :: eFrame
1343 +    real( kind = dp ), dimension(3,nLocal) :: f
1344 +    real( kind = dp ), dimension(3), intent(inout) :: fpair
1345 +
1346 +    real (kind = dp), dimension(3) :: ul1
1347 +    real (kind = dp), dimension(3) :: ul2
1348 +    real (kind = dp) :: dtdr
1349 +    real (kind = dp) :: dudx, dudy, dudz, u1dotu2
1350 +    integer :: me1, me2, id1, id2
1351 +    real (kind = dp) :: mu1, mu2
1352 +
1353 +    integer :: localError
1354 +
1355 +    if (.not.preRFCalculated) then
1356 +       call setReactionFieldPrefactor()
1357 +    endif
1358 +
1359 +    if (rij.le.rrf) then
1360 +
1361 +       if (rij.lt.rt) then
1362 +          dtdr = 0.0d0
1363 +       else
1364 +          !         write(*,*) 'rf correct in taper region'
1365 +          dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3)
1366 +       endif
1367 +
1368 + #ifdef IS_MPI
1369 +       me1 = atid_Row(atom1)
1370 +       ul1(1) = eFrame_Row(3,atom1)
1371 +       ul1(2) = eFrame_Row(6,atom1)
1372 +       ul1(3) = eFrame_Row(9,atom1)
1373 +
1374 +       me2 = atid_Col(atom2)
1375 +       ul2(1) = eFrame_Col(3,atom2)
1376 +       ul2(2) = eFrame_Col(6,atom2)
1377 +       ul2(3) = eFrame_Col(9,atom2)
1378 + #else
1379 +       me1 = atid(atom1)
1380 +       ul1(1) = eFrame(3,atom1)
1381 +       ul1(2) = eFrame(6,atom1)
1382 +       ul1(3) = eFrame(9,atom1)
1383 +
1384 +       me2 = atid(atom2)
1385 +       ul2(1) = eFrame(3,atom2)
1386 +       ul2(2) = eFrame(6,atom2)
1387 +       ul2(3) = eFrame(9,atom2)
1388 + #endif
1389 +
1390 +       mu1 = getDipoleMoment(me1)
1391 +       mu2 = getDipoleMoment(me2)
1392 +
1393 +       u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3)
1394 +
1395 +       dudx = - preRF*mu1*mu2*u1dotu2*dtdr*d(1)/rij
1396 +       dudy = - preRF*mu1*mu2*u1dotu2*dtdr*d(2)/rij
1397 +       dudz = - preRF*mu1*mu2*u1dotu2*dtdr*d(3)/rij
1398 +
1399 + #ifdef IS_MPI
1400 +       f_Row(1,atom1) = f_Row(1,atom1) + dudx
1401 +       f_Row(2,atom1) = f_Row(2,atom1) + dudy
1402 +       f_Row(3,atom1) = f_Row(3,atom1) + dudz
1403 +
1404 +       f_Col(1,atom2) = f_Col(1,atom2) - dudx
1405 +       f_Col(2,atom2) = f_Col(2,atom2) - dudy
1406 +       f_Col(3,atom2) = f_Col(3,atom2) - dudz
1407 + #else
1408 +       f(1,atom1) = f(1,atom1) + dudx
1409 +       f(2,atom1) = f(2,atom1) + dudy
1410 +       f(3,atom1) = f(3,atom1) + dudz
1411 +
1412 +       f(1,atom2) = f(1,atom2) - dudx
1413 +       f(2,atom2) = f(2,atom2) - dudy
1414 +       f(3,atom2) = f(3,atom2) - dudz
1415 + #endif
1416 +
1417 + #ifdef IS_MPI
1418 +       id1 = AtomRowToGlobal(atom1)
1419 +       id2 = AtomColToGlobal(atom2)
1420 + #else
1421 +       id1 = atom1
1422 +       id2 = atom2
1423 + #endif
1424 +
1425 +       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1426 +
1427 +          fpair(1) = fpair(1) + dudx
1428 +          fpair(2) = fpair(2) + dudy
1429 +          fpair(3) = fpair(3) + dudz
1430 +
1431 +       endif
1432 +
1433 +    end if
1434 +    return
1435 +  end subroutine rf_correct_forces
1436 +
1437   end module electrostatic_module

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines