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 2390 by chrisfen, Wed Oct 19 19:24:40 2005 UTC vs.
Revision 2394 by chrisfen, Sun Oct 23 21:08:08 2005 UTC

# Line 97 | Line 97 | module electrostatic_module
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 +  real(kind=dp), save :: preRF2 = 0.0_DP
101    logical, save :: preRFCalculated = .false.
102  
103   #ifdef __IFC
# Line 117 | Line 118 | module electrostatic_module
118    public :: doElectrostaticPair
119    public :: getCharge
120    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
122 >  public :: rf_self_self
123  
124    type :: Electrostatic
125       integer :: c_ident
# Line 176 | Line 173 | contains
173    subroutine setReactionFieldPrefactor
174      if (haveDefaultCutoff .and. haveDielectric) then
175         defaultCutoff2 = defaultCutoff*defaultCutoff
176 <       preRF = pre22 * 2.0d0*(dielectric-1.0d0) / &
176 >       preRF = (dielectric-1.0d0) / &
177              ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
178 +       preRF2 = 2.0d0*preRF
179         preRFCalculated = .true.
180      else if (.not.haveDefaultCutoff) then
181         call handleError("setReactionFieldPrefactor", "Default cutoff not set")
# Line 479 | Line 477 | contains
477      real (kind=dp) :: scale, sc2, bigR
478      real (kind=dp) :: varERFC, varEXP
479      real (kind=dp) :: limScale
480 +    real (kind=dp) :: preVal, rfVal
481  
482      if (.not.allocated(ElectrostaticMap)) then
483         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 487 | Line 486 | contains
486  
487      if (.not.summationMethodChecked) then
488         call checkSummationMethod()
490 –      
489      endif
490  
491 +    if (.not.preRFCalculated) then
492 +       call setReactionFieldPrefactor()
493 +    endif
494  
495   #ifdef IS_MPI
496      me1 = atid_Row(atom1)
# Line 643 | Line 644 | contains
644         if (j_is_Charge) then
645  
646            if (summationMethod .eq. UNDAMPED_WOLF) then
646 –
647               vterm = pre11 * q_i * q_j * (riji - rcuti)
648               vpair = vpair + vterm
649               epot = epot + sw*vterm
# Line 655 | Line 655 | contains
655               dudz = dudz + dudr * d(3)
656  
657            elseif (summationMethod .eq. DAMPED_WOLF) then
658 –
658               varERFC = derfc(dampingAlpha*rij)
659               varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
660               vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
# Line 671 | Line 670 | contains
670               dudy = dudy + dudr * d(2)
671               dudz = dudz + dudr * d(3)
672  
673 <          else
673 >          elseif (summationMethod .eq. REACTION_FIELD) then
674 >             preVal = pre11 * q_i * q_j
675 >             rfVal = preRF*rij*rij
676 >             vterm = preVal * ( riji + rfVal )
677 >             vpair = vpair + vterm
678 >             epot = epot + sw*vterm
679 >            
680 >             dudr  = sw * preVal * ( 2.0d0*rfVal - riji )*riji
681 >            
682 >             dudx = dudx + dudr * xhat
683 >             dudy = dudy + dudr * yhat
684 >             dudz = dudz + dudr * zhat
685  
686 +          else
687               vterm = pre11 * q_i * q_j * riji
688               vpair = vpair + vterm
689               epot = epot + sw*vterm
# Line 949 | Line 960 | contains
960                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
961               duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
962                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
963 +
964 +         elseif (summationMethod .eq. REACTION_FIELD) then
965 +             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
966 +
967 +             ri2 = riji * riji
968 +             ri3 = ri2 * riji
969 +             ri4 = ri2 * ri2
970 +
971 +             pref = pre22 * mu_i * mu_j
972 +              
973 +             vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
974 +                  preRF2*ct_ij )
975 +             vpair = vpair + vterm
976 +             epot = epot + sw*vterm
977 +            
978 +             a1 = 5.0d0 * ct_i * ct_j - ct_ij
979 +            
980 +             dudx = dudx + sw*pref*3.0d0*ri4 &
981 +                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
982 +             dudy = dudy + sw*pref*3.0d0*ri4 &
983 +                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
984 +             dudz = dudz + sw*pref*3.0d0*ri4 &
985 +                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
986 +            
987 +             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
988 +                  - preRF2*uz_j(1))
989 +             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
990 +                  - preRF2*uz_j(2))
991 +             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
992 +                  - preRF2*uz_j(3))
993 +             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
994 +                  - preRF2*uz_i(1))
995 +             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
996 +                  - preRF2*uz_i(2))
997 +             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
998 +                  - preRF2*uz_i(3))
999  
1000            else
1001               if (i_is_SplitDipole) then
# Line 1223 | Line 1270 | contains
1270  
1271    end subroutine destroyElectrostaticTypes
1272  
1273 <  subroutine accumulate_rf(atom1, atom2, rij, eFrame, taper)
1274 <
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 <
1273 >  subroutine rf_self_self(atom1, eFrame, rfpot, t, do_pot)
1274 >    logical, intent(in) :: do_pot
1275      integer, intent(in) :: atom1
1276 <    real(kind=dp), intent(in) :: mu1
1276 >    integer :: atid1
1277      real(kind=dp), dimension(9,nLocal) :: eFrame
1278 +    real(kind=dp), dimension(3,nLocal) :: t
1279 +    real(kind=dp) :: mu1
1280 +    real(kind=dp) :: preVal, epot, rfpot
1281 +    real(kind=dp) :: eix, eiy, eiz
1282  
1283 <    !! should work for both MPI and non-MPI version since this is not pairwise.
1284 <    rf(1,atom1) = rf(1,atom1) + eFrame(3,atom1)*mu1
1285 <    rf(2,atom1) = rf(2,atom1) + eFrame(6,atom1)*mu1
1286 <    rf(3,atom1) = rf(3,atom1) + eFrame(9,atom1)*mu1
1283 >    ! this is a local only array, so we use the local atom type id's:
1284 >    atid1 = atid(atom1)
1285 >    
1286 >    if (ElectrostaticMap(atid1)%is_Dipole) then
1287 >       mu1 = getDipoleMoment(atid1)
1288 >      
1289 >       preVal = pre22 * preRF2 * mu1*mu1
1290 >       rfpot = rfpot - 0.5d0*preVal
1291  
1292 <    return
1293 <  end subroutine accumulate_self_rf
1292 >       ! The self-correction term adds into the reaction field vector
1293 >      
1294 >       eix = preVal * eFrame(3,atom1)
1295 >       eiy = preVal * eFrame(6,atom1)
1296 >       eiz = preVal * eFrame(9,atom1)
1297  
1298 <  subroutine reaction_field_final(a1, mu1, eFrame, rfpot, t, do_pot)
1298 >       ! once again, this is self-self, so only the local arrays are needed
1299 >       ! even for MPI jobs:
1300 >      
1301 >       t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1302 >            eFrame(9,atom1)*eiy
1303 >       t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1304 >            eFrame(3,atom1)*eiz
1305 >       t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1306 >            eFrame(6,atom1)*eix
1307  
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()
1308      endif
1309 <
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 <
1309 >    
1310      return
1311 <  end subroutine reaction_field_final
1311 >  end subroutine rf_self_self
1312  
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 –
1313   end module electrostatic_module

Diff Legend

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