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 2095 by gezelter, Wed Mar 9 15:44:59 2005 UTC vs.
Revision 2118 by gezelter, Fri Mar 11 15:53:18 2005 UTC

# Line 54 | Line 54 | module electrostatic_module
54  
55    PRIVATE
56  
57 +  !! these prefactors convert the multipole interactions into kcal / mol
58 +  !! all were computed assuming distances are measured in angstroms
59 +  !! Charge-Charge, assuming charges are measured in electrons
60    real(kind=dp), parameter :: pre11 = 332.0637778_dp
61 <  real(kind=dp), parameter :: pre12 = 69.13291783_dp
62 <  real(kind=dp), parameter :: pre22 = 14.39289874_dp
61 >  !! Charge-Dipole, assuming charges are measured in electrons, and
62 >  !! dipoles are measured in debyes
63 >  real(kind=dp), parameter :: pre12 = 69.13373_dp
64 >  !! Dipole-Dipole, assuming dipoles are measured in debyes
65 >  real(kind=dp), parameter :: pre22 = 14.39325_dp
66 >  !! Charge-Quadrupole, assuming charges are measured in electrons, and
67 >  !! quadrupoles are measured in 10^-26 esu cm^2
68 >  !! This unit is also known affectionately as an esu centi-barn.
69 >  real(kind=dp), parameter :: pre14 = 69.13373_dp
70  
71    public :: newElectrostaticType
72    public :: setCharge
# Line 317 | Line 327 | contains
327      integer :: me1, me2, id1, id2
328      real (kind=dp) :: q_i, q_j, mu_i, mu_j, d_i, d_j
329      real (kind=dp) :: ct_i, ct_j, ct_ij, a1
330 <    real (kind=dp) :: riji, ri2, ri3, ri4
330 >    real (kind=dp) :: riji, ri, ri2, ri3, ri4
331      real (kind=dp) :: pref, vterm, epot, dudr    
332 +    real (kind=dp) :: xhat, yhat, zhat
333      real (kind=dp) :: dudx, dudy, dudz
334      real (kind=dp) :: drdxj, drdyj, drdzj
335      real (kind=dp) :: duduix, duduiy, duduiz, dudujx, dudujy, dudujz
336 +    real (kind=dp) :: scale, sc2, bigR
337  
326 –
338      if (.not.allocated(ElectrostaticMap)) then
339         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
340         return
# Line 341 | Line 352 | contains
352  
353      riji = 1.0d0 / rij
354  
355 <    !! these are also useful as the unit vector of \vec{r}
356 <    !! \hat{r} = \vec{r} / r =   {(x_j-x_i) / r, (y_j-y_i)/r, (z_j-z_i)/r}
355 >    xhat = d(1) * riji
356 >    yhat = d(2) * riji
357 >    zhat = d(3) * riji
358  
359 <    drdxj = d(1) * riji
360 <    drdyj = d(2) * riji
361 <    drdzj = d(3) * riji
359 >    drdxj = xhat
360 >    drdyj = yhat
361 >    drdzj = zhat
362  
363      !! logicals
364  
# Line 436 | Line 448 | contains
448  
449         if (j_is_Dipole) then
450  
451 <          ri2 = riji * riji
452 <          ri3 = ri2 * riji
451 >          if (j_is_SplitDipole) then
452 >             BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
453 >             ri = 1.0_dp / BigR
454 >             scale = rij * ri
455 >          else
456 >             ri = riji
457 >             scale = 1.0_dp
458 >          endif
459  
460 +          ri2 = ri * ri
461 +          ri3 = ri2 * ri
462 +          sc2 = scale * scale
463 +            
464            pref = pre12 * q_i * mu_j
465 <          vterm = pref * ct_j * riji * riji
465 >          vterm = pref * ct_j * ri2 * scale
466            vpair = vpair + vterm
467            epot = epot + sw * vterm
468  
469 <          dudx = dudx + pref * sw * ri3 * ( ul_j(1) + 3.0d0 * ct_j * drdxj)
470 <          dudy = dudy + pref * sw * ri3 * ( ul_j(2) + 3.0d0 * ct_j * drdyj)
471 <          dudz = dudz + pref * sw * ri3 * ( ul_j(3) + 3.0d0 * ct_j * drdzj)
469 >          !! this has a + sign in the () because the rij vector is
470 >          !! r_j - r_i and the charge-dipole potential takes the origin
471 >          !! as the point dipole, which is atom j in this case.
472  
473 <          dudujx = dudujx - pref * sw * ri2 * drdxj
474 <          dudujy = dudujy - pref * sw * ri2 * drdyj
475 <          dudujz = dudujz - pref * sw * ri2 * drdzj
473 >          dudx = dudx + pref * sw * ri3 * ( ul_j(1) + 3.0d0*ct_j*xhat*sc2)
474 >          dudy = dudy + pref * sw * ri3 * ( ul_j(2) + 3.0d0*ct_j*yhat*sc2)
475 >          dudz = dudz + pref * sw * ri3 * ( ul_j(3) + 3.0d0*ct_j*zhat*sc2)
476 >
477 >          dudujx = dudujx - pref * sw * ri2 * xhat * scale
478 >          dudujy = dudujy - pref * sw * ri2 * yhat * scale
479 >          dudujz = dudujz - pref * sw * ri2 * zhat * scale
480            
481         endif
482 +
483      endif
484    
485      if (i_is_Dipole) then
486        
487         if (j_is_Charge) then
488  
489 <          ri2 = riji * riji
490 <          ri3 = ri2 * riji
489 >          if (i_is_SplitDipole) then
490 >             BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
491 >             ri = 1.0_dp / BigR
492 >             scale = rij * ri
493 >          else
494 >             ri = riji
495 >             scale = 1.0_dp
496 >          endif
497  
498 +          ri2 = ri * ri
499 +          ri3 = ri2 * ri
500 +          sc2 = scale * scale
501 +            
502            pref = pre12 * q_j * mu_i
503 <          vterm = pref * ct_i * riji * riji
503 >          vterm = pref * ct_i * ri2 * scale
504            vpair = vpair + vterm
505            epot = epot + sw * vterm
506  
507 <          dudx = dudx + pref * sw * ri3 * ( ul_i(1) - 3.0d0 * ct_i * drdxj)
508 <          dudy = dudy + pref * sw * ri3 * ( ul_i(2) - 3.0d0 * ct_i * drdyj)
509 <          dudz = dudz + pref * sw * ri3 * ( ul_i(3) - 3.0d0 * ct_i * drdzj)
507 >          dudx = dudx + pref * sw * ri3 * ( ul_i(1) - 3.0d0 * ct_i * xhat*sc2)
508 >          dudy = dudy + pref * sw * ri3 * ( ul_i(2) - 3.0d0 * ct_i * yhat*sc2)
509 >          dudz = dudz + pref * sw * ri3 * ( ul_i(3) - 3.0d0 * ct_i * zhat*sc2)
510  
511 <          duduix = duduix + pref * sw * ri2 * drdxj
512 <          duduiy = duduiy + pref * sw * ri2 * drdyj
513 <          duduiz = duduiz + pref * sw * ri2 * drdzj
511 >          duduix = duduix + pref * sw * ri2 * xhat * scale
512 >          duduiy = duduiy + pref * sw * ri2 * yhat * scale
513 >          duduiz = duduiz + pref * sw * ri2 * zhat * scale
514         endif
515  
516         if (j_is_Dipole) then
517  
518 +          if (i_is_SplitDipole) then
519 +             if (j_is_SplitDipole) then
520 +                BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
521 +             else
522 +                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
523 +             endif
524 +             ri = 1.0_dp / BigR
525 +             scale = rij * ri                
526 +          else
527 +             if (j_is_SplitDipole) then
528 +                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
529 +                ri = 1.0_dp / BigR
530 +                scale = rij * ri                            
531 +             else                
532 +                ri = riji
533 +                scale = 1.0_dp
534 +             endif
535 +          endif
536 +
537            ct_ij = ul_i(1)*ul_j(1) + ul_i(2)*ul_j(2) + ul_i(3)*ul_j(3)
538 <          ri2 = riji * riji
539 <          ri3 = ri2 * riji
538 >
539 >          ri2 = ri * ri
540 >          ri3 = ri2 * ri
541            ri4 = ri2 * ri2
542 +          sc2 = scale * scale
543  
544            pref = pre22 * mu_i * mu_j
545 <          vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j)
545 >          vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j * sc2)
546            vpair = vpair + vterm
547            epot = epot + sw * vterm
548            
549 <          a1 = 5.0d0 * ct_i * ct_j - ct_ij
549 >          a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
550  
551 <          dudx = dudx + pref*sw*3.0d0*ri4*(a1*drdxj-ct_i*ul_j(1)-ct_j*ul_i(1))
552 <          dudy = dudy + pref*sw*3.0d0*ri4*(a1*drdyj-ct_i*ul_j(2)-ct_j*ul_i(2))
553 <          dudz = dudz + pref*sw*3.0d0*ri4*(a1*drdzj-ct_i*ul_j(3)-ct_j*ul_i(3))
551 >          dudx=dudx+pref*sw*3.0d0*ri4*scale*(a1*xhat-ct_i*ul_j(1)-ct_j*ul_i(1))
552 >          dudy=dudy+pref*sw*3.0d0*ri4*scale*(a1*yhat-ct_i*ul_j(2)-ct_j*ul_i(2))
553 >          dudz=dudz+pref*sw*3.0d0*ri4*scale*(a1*zhat-ct_i*ul_j(3)-ct_j*ul_i(3))
554  
555 <          duduix = duduix + pref*sw*ri3*(ul_j(1) - 3.0d0*ct_j*drdxj)
556 <          duduiy = duduiy + pref*sw*ri3*(ul_j(2) - 3.0d0*ct_j*drdyj)
557 <          duduiz = duduiz + pref*sw*ri3*(ul_j(3) - 3.0d0*ct_j*drdzj)
555 >          duduix = duduix + pref*sw*ri3*(ul_j(1) - 3.0d0*ct_j*xhat*sc2)
556 >          duduiy = duduiy + pref*sw*ri3*(ul_j(2) - 3.0d0*ct_j*yhat*sc2)
557 >          duduiz = duduiz + pref*sw*ri3*(ul_j(3) - 3.0d0*ct_j*zhat*sc2)
558  
559 <          dudujx = dudujx + pref*sw*ri3*(ul_i(1) - 3.0d0*ct_i*drdxj)
560 <          dudujy = dudujy + pref*sw*ri3*(ul_i(2) - 3.0d0*ct_i*drdyj)
561 <          dudujz = dudujz + pref*sw*ri3*(ul_i(3) - 3.0d0*ct_i*drdzj)
559 >          dudujx = dudujx + pref*sw*ri3*(ul_i(1) - 3.0d0*ct_i*xhat*sc2)
560 >          dudujy = dudujy + pref*sw*ri3*(ul_i(2) - 3.0d0*ct_i*yhat*sc2)
561 >          dudujz = dudujz + pref*sw*ri3*(ul_i(3) - 3.0d0*ct_i*zhat*sc2)
562         endif
563  
564      endif
# Line 577 | Line 635 | end module electrostatic_module
635    end subroutine doElectrostaticPair
636    
637   end module electrostatic_module
580 –

Diff Legend

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