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 2398 by chrisfen, Mon Oct 24 14:06:36 2005 UTC vs.
Revision 2399 by chrisfen, Wed Oct 26 23:31:40 2005 UTC

# Line 120 | Line 120 | module electrostatic_module
120    public :: getDipoleMoment
121    public :: destroyElectrostaticTypes
122    public :: rf_self_self
123 +  public :: rf_self_excludes
124  
125    type :: Electrostatic
126       integer :: c_ident
# Line 437 | Line 438 | contains
438  
439  
440    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
441 <       vpair, fpair, pot, eFrame, f, t, do_pot)
441 >       vpair, fpair, pot, eFrame, f, t, do_pot, indirect_only)
442  
443 <    logical, intent(in) :: do_pot
443 >    logical, intent(in) :: do_pot, indirect_only
444  
445      integer, intent(in) :: atom1, atom2
446      integer :: localError
# Line 674 | Line 675 | contains
675               preVal = pre11 * q_i * q_j
676               rfVal = preRF*rij*rij
677               vterm = preVal * ( riji + rfVal )
678 +            
679               vpair = vpair + vterm
680               epot = epot + sw*vterm
681              
# Line 727 | Line 729 | contains
729               duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
730  
731            elseif (summationMethod .eq. REACTION_FIELD) then
732 <             ri2 = ri * ri
733 <             ri3 = ri2 * ri
732 >             ri2 = riji * riji
733 >             ri3 = ri2 * riji
734      
735               pref = pre12 * q_i * mu_j
736               vterm = - pref * ct_j * ( ri2 - preRF2*rij )
# Line 899 | Line 901 | contains
901               vterm = pref * ct_i * (ri2 - rcuti2)
902               vpair = vpair + vterm
903               epot = epot + sw*vterm
902 –            
903 –             !! this has a + sign in the () because the rij vector is
904 –             !! r_j - r_i and the charge-dipole potential takes the origin
905 –             !! as the point dipole, which is atom j in this case.
904              
905               dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
906                    - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
# Line 911 | Line 909 | contains
909               dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
910                    - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
911              
912 <             duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
913 <             duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
914 <             duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
912 >             duduz_i(1) = duduz_i(1) + sw*pref*( ri2*xhat - d(1)*rcuti3 )
913 >             duduz_i(2) = duduz_i(2) + sw*pref*( ri2*yhat - d(2)*rcuti3 )
914 >             duduz_i(3) = duduz_i(3) + sw*pref*( ri2*zhat - d(3)*rcuti3 )
915  
916            elseif (summationMethod .eq. REACTION_FIELD) then
917 <             ri2 = ri * ri
918 <             ri3 = ri2 * ri
917 >             ri2 = riji * riji
918 >             ri3 = ri2 * riji
919  
920               pref = pre12 * q_j * mu_i
921 <             vterm = pref * ct_i * ( ri2 - preRF*rij )
921 >             vterm = pref * ct_i * ( ri2 - preRF2*rij )
922               vpair = vpair + vterm
923               epot = epot + sw*vterm
924              
925 <             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0*ct_i*xhat - &
926 <                                             preRF*uz_i(1) )
927 <             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0*ct_i*yhat - &
928 <                                             preRF*uz_i(2) )
929 <             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0*ct_i*zhat - &
930 <                                             preRF*uz_i(3) )
925 >             dudx = dudx + sw*pref * ( ri3*(uz_i(1) - 3.0d0*ct_i*xhat) - &
926 >                  preRF2*uz_i(1) )
927 >             dudy = dudy + sw*pref * ( ri3*(uz_i(2) - 3.0d0*ct_i*yhat) - &
928 >                  preRF2*uz_i(2) )
929 >             dudz = dudz + sw*pref * ( ri3*(uz_i(3) - 3.0d0*ct_i*zhat) - &
930 >                  preRF2*uz_i(3) )
931              
932 <             duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF*rij )
933 <             duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF*rij )
934 <             duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF*rij )
932 >             duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF2*rij )
933 >             duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF2*rij )
934 >             duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF2*rij )
935  
936            else
937               if (i_is_SplitDipole) then
# Line 1322 | Line 1320 | contains
1320      real(kind=dp) :: mu1
1321      real(kind=dp) :: preVal, epot, rfpot
1322      real(kind=dp) :: eix, eiy, eiz
1323 +
1324 +    if (.not.preRFCalculated) then
1325 +       call setReactionFieldPrefactor()
1326 +    endif
1327  
1328      ! this is a local only array, so we use the local atom type id's:
1329      atid1 = atid(atom1)
# Line 1331 | Line 1333 | contains
1333        
1334         preVal = pre22 * preRF2 * mu1*mu1
1335         rfpot = rfpot - 0.5d0*preVal
1336 <
1336 >      
1337         ! The self-correction term adds into the reaction field vector
1338        
1339         eix = preVal * eFrame(3,atom1)
1340         eiy = preVal * eFrame(6,atom1)
1341         eiz = preVal * eFrame(9,atom1)
1342 <
1342 >      
1343         ! once again, this is self-self, so only the local arrays are needed
1344         ! even for MPI jobs:
1345        
# Line 1347 | Line 1349 | contains
1349              eFrame(3,atom1)*eiz
1350         t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1351              eFrame(6,atom1)*eix
1352 <
1352 >      
1353      endif
1354      
1355      return
1356    end subroutine rf_self_self
1357  
1358 +  subroutine rf_self_excludes(atom1, atom2, sw, eFrame, d, rij, vpair, rfpot, &
1359 +       f, t, do_pot)
1360 +    logical, intent(in) :: do_pot
1361 +    integer, intent(in) :: atom1
1362 +    integer, intent(in) :: atom2
1363 +    logical :: i_is_Charge, j_is_Charge
1364 +    logical :: i_is_Dipole, j_is_Dipole
1365 +    integer :: atid1
1366 +    integer :: atid2
1367 +    real(kind=dp), intent(in) :: rij
1368 +    real(kind=dp), intent(in) :: sw
1369 +    real(kind=dp), intent(in), dimension(3) :: d
1370 +    real(kind=dp), intent(inout) :: vpair
1371 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1372 +    real(kind=dp), dimension(3,nLocal) :: f
1373 +    real(kind=dp), dimension(3,nLocal) :: t
1374 +    real (kind = dp), dimension(3) :: duduz_i
1375 +    real (kind = dp), dimension(3) :: duduz_j
1376 +    real (kind = dp), dimension(3) :: uz_i
1377 +    real (kind = dp), dimension(3) :: uz_j
1378 +    real(kind=dp) :: q_i, q_j, mu_i, mu_j
1379 +    real(kind=dp) :: xhat, yhat, zhat
1380 +    real(kind=dp) :: ct_i, ct_j
1381 +    real(kind=dp) :: ri2, ri3, riji, vterm
1382 +    real(kind=dp) :: pref, preVal, rfVal, rfpot
1383 +    real(kind=dp) :: dudx, dudy, dudz, dudr
1384 +
1385 +    if (.not.preRFCalculated) then
1386 +       call setReactionFieldPrefactor()
1387 +    endif
1388 +
1389 +    dudx = 0.0d0
1390 +    dudy = 0.0d0
1391 +    dudz = 0.0d0
1392 +
1393 +    riji = 1.0d0/rij
1394 +
1395 +    xhat = d(1) * riji
1396 +    yhat = d(2) * riji
1397 +    zhat = d(3) * riji
1398 +
1399 +    ! this is a local only array, so we use the local atom type id's:
1400 +    atid1 = atid(atom1)
1401 +    atid2 = atid(atom2)
1402 +    i_is_Charge = ElectrostaticMap(atid1)%is_Charge
1403 +    j_is_Charge = ElectrostaticMap(atid2)%is_Charge
1404 +    i_is_Dipole = ElectrostaticMap(atid1)%is_Dipole
1405 +    j_is_Dipole = ElectrostaticMap(atid2)%is_Dipole
1406 +
1407 +    if (i_is_Charge.and.j_is_Charge) then
1408 +       q_i = ElectrostaticMap(atid1)%charge
1409 +       q_j = ElectrostaticMap(atid2)%charge
1410 +      
1411 +       preVal = pre11 * q_i * q_j
1412 +       rfVal = preRF*rij*rij
1413 +       vterm = preVal * rfVal
1414 +
1415 +       rfpot = rfpot + sw*vterm
1416 +      
1417 +       dudr  = sw*preVal * 2.0d0*rfVal*riji
1418 +              
1419 +       dudx = dudx + dudr * xhat
1420 +       dudy = dudy + dudr * yhat
1421 +       dudz = dudz + dudr * zhat
1422 +
1423 +    elseif (i_is_Charge.and.j_is_Dipole) then
1424 +       q_i = ElectrostaticMap(atid1)%charge
1425 +       mu_j = ElectrostaticMap(atid2)%dipole_moment
1426 +       uz_j(1) = eFrame(3,atom2)
1427 +       uz_j(2) = eFrame(6,atom2)
1428 +       uz_j(3) = eFrame(9,atom2)
1429 +       ct_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
1430 +
1431 +       ri2 = riji * riji
1432 +       ri3 = ri2 * riji
1433 +      
1434 +       pref = pre12 * q_i * mu_j
1435 +       vterm = - pref * ct_j * ( ri2 - preRF2*rij )
1436 +       rfpot = rfpot + sw*vterm
1437 +
1438 +       dudx = dudx - sw*pref*( ri3*(uz_j(1)-3.0d0*ct_j*xhat) - preRF2*uz_j(1) )
1439 +       dudy = dudy - sw*pref*( ri3*(uz_j(2)-3.0d0*ct_j*yhat) - preRF2*uz_j(2) )
1440 +       dudz = dudz - sw*pref*( ri3*(uz_j(3)-3.0d0*ct_j*zhat) - preRF2*uz_j(3) )
1441 +
1442 +       duduz_j(1) = duduz_j(1) - sw * pref * xhat * ( ri2 - preRF2*rij )
1443 +       duduz_j(2) = duduz_j(2) - sw * pref * yhat * ( ri2 - preRF2*rij )
1444 +       duduz_j(3) = duduz_j(3) - sw * pref * zhat * ( ri2 - preRF2*rij )
1445 +              
1446 +    elseif (i_is_Dipole.and.j_is_Charge) then
1447 +       mu_i = ElectrostaticMap(atid1)%dipole_moment
1448 +       q_j = ElectrostaticMap(atid2)%charge
1449 +       uz_i(1) = eFrame(3,atom1)
1450 +       uz_i(2) = eFrame(6,atom1)
1451 +       uz_i(3) = eFrame(9,atom1)
1452 +       ct_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
1453 +
1454 +       ri2 = riji * riji
1455 +       ri3 = ri2 * riji
1456 +      
1457 +       pref = pre12 * q_j * mu_i
1458 +       vterm = pref * ct_i * ( ri2 - preRF2*rij )
1459 +       rfpot = rfpot + sw*vterm
1460 +      
1461 +       dudx = dudx + sw*pref*( ri3*(uz_i(1)-3.0d0*ct_i*xhat) - preRF2*uz_i(1) )
1462 +       dudy = dudy + sw*pref*( ri3*(uz_i(2)-3.0d0*ct_i*yhat) - preRF2*uz_i(2) )
1463 +       dudz = dudz + sw*pref*( ri3*(uz_i(3)-3.0d0*ct_i*zhat) - preRF2*uz_i(3) )
1464 +      
1465 +       duduz_i(1) = duduz_i(1) + sw * pref * xhat * ( ri2 - preRF2*rij )
1466 +       duduz_i(2) = duduz_i(2) + sw * pref * yhat * ( ri2 - preRF2*rij )
1467 +       duduz_i(3) = duduz_i(3) + sw * pref * zhat * ( ri2 - preRF2*rij )
1468 +      
1469 +    endif
1470 +    
1471 +    ! accumulate the forces and torques resulting from the RF self term
1472 +    f(1,atom1) = f(1,atom1) + dudx
1473 +    f(2,atom1) = f(2,atom1) + dudy
1474 +    f(3,atom1) = f(3,atom1) + dudz
1475 +    
1476 +    f(1,atom2) = f(1,atom2) - dudx
1477 +    f(2,atom2) = f(2,atom2) - dudy
1478 +    f(3,atom2) = f(3,atom2) - dudz
1479 +    
1480 +    if (i_is_Dipole) then
1481 +       t(1,atom1)=t(1,atom1) - uz_i(2)*duduz_i(3) + uz_i(3)*duduz_i(2)
1482 +       t(2,atom1)=t(2,atom1) - uz_i(3)*duduz_i(1) + uz_i(1)*duduz_i(3)
1483 +       t(3,atom1)=t(3,atom1) - uz_i(1)*duduz_i(2) + uz_i(2)*duduz_i(1)
1484 +    elseif (j_is_Dipole) then
1485 +       t(1,atom2)=t(1,atom2) - uz_j(2)*duduz_j(3) + uz_j(3)*duduz_j(2)
1486 +       t(2,atom2)=t(2,atom2) - uz_j(3)*duduz_j(1) + uz_j(1)*duduz_j(3)
1487 +       t(3,atom2)=t(3,atom2) - uz_j(1)*duduz_j(2) + uz_j(2)*duduz_j(1)
1488 +    endif
1489 +
1490 +    return
1491 +  end subroutine rf_self_excludes
1492 +
1493   end module electrostatic_module

Diff Legend

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