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 2402 by chrisfen, Tue Nov 1 19:09:30 2005 UTC vs.
Revision 2411 by chrisfen, Wed Nov 2 21:01:21 2005 UTC

# Line 58 | Line 58 | module electrostatic_module
58   #define __FORTRAN90
59   #include "UseTheForce/DarkSide/fInteractionMap.h"
60   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61 + #include "UseTheForce/DarkSide/fScreeningMethod.h"
62  
63  
64    !! these prefactors convert the multipole interactions into kcal / mol
# Line 74 | Line 75 | module electrostatic_module
75    !! This unit is also known affectionately as an esu centi-barn.
76    real(kind=dp), parameter :: pre14 = 69.13373_dp
77  
78 <  !! variables to handle different summation methods for long-range electrostatics:
78 >  !! variables to handle different summation methods for long-range
79 >  !! electrostatics:
80    integer, save :: summationMethod = NONE
81 +  integer, save :: screeningMethod = UNDAMPED
82    logical, save :: summationMethodChecked = .false.
83    real(kind=DP), save :: defaultCutoff = 0.0_DP
84    real(kind=DP), save :: defaultCutoff2 = 0.0_DP
# Line 104 | Line 107 | module electrostatic_module
107   #endif
108    
109    public :: setElectrostaticSummationMethod
110 +  public :: setScreeningMethod
111    public :: setElectrostaticCutoffRadius
112 <  public :: setDampedWolfAlpha
112 >  public :: setDampingAlpha
113    public :: setReactionFieldDielectric
114    public :: newElectrostaticType
115    public :: setCharge
# Line 147 | Line 151 | contains
151  
152    end subroutine setElectrostaticSummationMethod
153  
154 +  subroutine setScreeningMethod(the_SM)
155 +    integer, intent(in) :: the_SM    
156 +    screeningMethod = the_SM
157 +  end subroutine setScreeningMethod
158 +
159    subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
160      real(kind=dp), intent(in) :: thisRcut
161      real(kind=dp), intent(in) :: thisRsw
# Line 156 | Line 165 | contains
165      haveDefaultCutoff = .true.
166    end subroutine setElectrostaticCutoffRadius
167  
168 <  subroutine setDampedWolfAlpha(thisAlpha)
168 >  subroutine setDampingAlpha(thisAlpha)
169      real(kind=dp), intent(in) :: thisAlpha
170      dampingAlpha = thisAlpha
171      haveDampingAlpha = .true.
172 <  end subroutine setDampedWolfAlpha
172 >  end subroutine setDampingAlpha
173    
174    subroutine setReactionFieldDielectric(thisDielectric)
175      real(kind=dp), intent(in) :: thisDielectric
# Line 389 | Line 398 | contains
398      rcuti3 = rcuti2*rcuti
399      rcuti4 = rcuti2*rcuti2
400  
401 <    if (summationMethod .eq. DAMPED_WOLF) then
401 >    if (screeningMethod .eq. DAMPED) then
402         if (.not.haveDampingAlpha) then
403            call handleError("checkSummationMethod", "no Damping Alpha set!")
404         endif
# Line 420 | Line 429 | contains
429      summationMethodChecked = .true.
430    end subroutine checkSummationMethod
431  
432 < !!$
424 < !!$  subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
425 < !!$       vpair, fpair, pot, eFrame, f, t, do_pot)
432 >
433    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
434 <       vpair, fpair, pot, eFrame, f, t, do_pot, felec)
434 >       vpair, fpair, pot, eFrame, f, t, do_pot)
435  
436      logical, intent(in) :: do_pot
437  
# Line 435 | Line 442 | contains
442      real(kind=dp), intent(in), dimension(3) :: d
443      real(kind=dp), intent(inout) :: vpair
444      real(kind=dp), intent(inout), dimension(3) :: fpair    
438 –    real(kind=dp), intent(inout), dimension(3) :: felec
445  
446      real( kind = dp ) :: pot
447      real( kind = dp ), dimension(9,nLocal) :: eFrame
448      real( kind = dp ), dimension(3,nLocal) :: f
449 +    real( kind = dp ), dimension(3,nLocal) :: felec
450      real( kind = dp ), dimension(3,nLocal) :: t
451  
452      real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
# Line 484 | Line 491 | contains
491      me2 = atid(atom2)
492   #endif
493  
494 + !!$    if (rij .ge. defaultCutoff) then
495 + !!$       write(*,*) 'warning: rij = ', rij, ' rcut = ', defaultCutoff, ' sw = ', sw
496 + !!$    endif
497 +
498      !! some variables we'll need independent of electrostatic type:
499  
500      riji = 1.0d0 / rij
# Line 627 | Line 638 | contains
638  
639         if (j_is_Charge) then
640  
641 <          if (summationMethod .eq. UNDAMPED_WOLF) then
641 >          if (summationMethod .eq. SHIFTED_POTENTIAL) then
642 >
643               vterm = pre11 * q_i * q_j * (riji - rcuti)
644               vpair = vpair + vterm
645               epot = epot + sw*vterm
646              
647 <             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)
647 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji)
648              
649               dudx = dudx + dudr * xhat
650               dudy = dudy + dudr * yhat
651               dudz = dudz + dudr * zhat
652  
653 <          elseif (summationMethod .eq. DAMPED_WOLF) then
654 <             varERFC = derfc(dampingAlpha*rij)
655 <             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
656 <             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
653 > !!$          elseif (summationMethod .eq. DAMPED_WOLF) then
654 > !!$             varERFC = derfc(dampingAlpha*rij)
655 > !!$             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
656 > !!$             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
657 > !!$             vpair = vpair + vterm
658 > !!$             epot = epot + sw*vterm
659 > !!$            
660 > !!$!             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
661 > !!$!                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
662 > !!$!                  + alphaPi*constEXP*rcuti)) )
663 > !!$             dudr  = -sw*pre11*q_i*q_j * (varERFC*riji*riji &
664 > !!$                  + alphaPi*varEXP*riji)
665 > !!$            
666 > !!$             dudx = dudx + dudr * xhat
667 > !!$             dudy = dudy + dudr * yhat
668 > !!$             dudz = dudz + dudr * zhat
669 >
670 >          elseif (summationMethod .eq. SHIFTED_FORCE) then
671 >             vterm = pre11 * q_i * q_j * (riji + rij*rcuti2 - 2.0d0*rcuti)
672               vpair = vpair + vterm
673               epot = epot + sw*vterm
674              
675 <             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
676 <                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
650 <                  + alphaPi*constEXP*rcuti)) )
651 <            
675 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)
676 >                          
677               dudx = dudx + dudr * xhat
678               dudy = dudy + dudr * yhat
679               dudz = dudz + dudr * zhat
680  
681 + !!$          elseif (summationMethod .eq. DAMPED_SHIFTED_FORCE) then
682 + !!$             varERFC = derfc(dampingAlpha*rij)
683 + !!$             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
684 + !!$             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
685 + !!$             vpair = vpair + vterm
686 + !!$             epot = epot + sw*vterm
687 + !!$            
688 + !!$!             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
689 + !!$!                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
690 + !!$!                  + alphaPi*constEXP*rcuti)) )
691 + !!$             dudr  = -sw*pre11*q_i*q_j * (varERFC*riji*riji &
692 + !!$                  + alphaPi*varEXP*riji)
693 + !!$            
694 + !!$             dudx = dudx + dudr * xhat
695 + !!$             dudy = dudy + dudr * yhat
696 + !!$             dudz = dudz + dudr * zhat
697 +
698            elseif (summationMethod .eq. REACTION_FIELD) then
699               preVal = pre11 * q_i * q_j
700               rfVal = preRF*rij*rij
# Line 686 | Line 728 | contains
728  
729            pref = pre12 * q_i * mu_j
730  
731 <          if (summationMethod .eq. UNDAMPED_WOLF) then
732 <             ri2 = riji * riji
733 <             ri3 = ri2 * riji
734 <
735 <             pref = pre12 * q_i * mu_j
736 <             vterm = - pref * ct_j * (ri2 - rcuti2)
737 <             vpair = vpair + vterm
738 <             epot = epot + sw*vterm
739 <            
740 <             !! this has a + sign in the () because the rij vector is
741 <             !! r_j - r_i and the charge-dipole potential takes the origin
742 <             !! as the point dipole, which is atom j in this case.
743 <            
744 <             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
745 <                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
746 <             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
747 <                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
748 <             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
749 <                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
750 <            
751 <             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
752 <             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
753 <             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
731 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
732 > !!$             ri2 = riji * riji
733 > !!$             ri3 = ri2 * riji
734 > !!$
735 > !!$             pref = pre12 * q_i * mu_j
736 > !!$             vterm = - pref * ct_j * (ri2 - rcuti2)
737 > !!$             vpair = vpair + vterm
738 > !!$             epot = epot + sw*vterm
739 > !!$            
740 > !!$             !! this has a + sign in the () because the rij vector is
741 > !!$             !! r_j - r_i and the charge-dipole potential takes the origin
742 > !!$             !! as the point dipole, which is atom j in this case.
743 > !!$            
744 > !!$             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
745 > !!$                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
746 > !!$             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
747 > !!$                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
748 > !!$             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
749 > !!$                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
750 > !!$            
751 > !!$             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
752 > !!$             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
753 > !!$             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
754 > !!$
755 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
756  
757 <          elseif (summationMethod .eq. REACTION_FIELD) then
757 >          if (summationMethod .eq. REACTION_FIELD) then
758               ri2 = riji * riji
759               ri3 = ri2 * riji
760      
# Line 775 | Line 819 | contains
819            cy2 = cy_j * cy_j
820            cz2 = cz_j * cz_j
821  
822 <          if (summationMethod .eq. UNDAMPED_WOLF) then
823 <             pref =  pre14 * q_i / 3.0_dp
824 <             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
825 <                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
826 <                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
827 <             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
828 <                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
829 <                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
830 <             vpair = vpair + ( vterm1 - vterm2 )
831 <             epot = epot + sw*( vterm1 - vterm2 )
832 <            
833 <             dudx = dudx - (5.0_dp * &
834 <                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
835 <                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
836 <                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
837 <                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
838 <                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
839 <                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
840 <                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
841 <             dudy = dudy - (5.0_dp * &
842 <                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
843 <                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
844 <                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
845 <                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
846 <                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
847 <                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
848 <                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
849 <             dudz = dudz - (5.0_dp * &
850 <                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
851 <                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
852 <                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
853 <                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
854 <                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
855 <                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
856 <                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
857 <            
858 <             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
859 <                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
860 <             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
861 <                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
862 <             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
863 <                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
864 <            
865 <             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
866 <                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
867 <             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
868 <                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
869 <             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
870 <                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
871 <            
872 <             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
873 <                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
874 <             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
875 <                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
876 <             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
877 <                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
878 <        
879 <          else
822 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
823 > !!$             pref =  pre14 * q_i / 3.0_dp
824 > !!$             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
825 > !!$                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
826 > !!$                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
827 > !!$             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
828 > !!$                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
829 > !!$                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
830 > !!$             vpair = vpair + ( vterm1 - vterm2 )
831 > !!$             epot = epot + sw*( vterm1 - vterm2 )
832 > !!$            
833 > !!$             dudx = dudx - (5.0_dp * &
834 > !!$                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
835 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
836 > !!$                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
837 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
838 > !!$                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
839 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
840 > !!$                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
841 > !!$             dudy = dudy - (5.0_dp * &
842 > !!$                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
843 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
844 > !!$                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
845 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
846 > !!$                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
847 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
848 > !!$                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
849 > !!$             dudz = dudz - (5.0_dp * &
850 > !!$                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
851 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
852 > !!$                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
853 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
854 > !!$                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
855 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
856 > !!$                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
857 > !!$            
858 > !!$             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
859 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
860 > !!$             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
861 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
862 > !!$             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
863 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
864 > !!$            
865 > !!$             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
866 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
867 > !!$             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
868 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
869 > !!$             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
870 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
871 > !!$            
872 > !!$             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
873 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
874 > !!$             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
875 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
876 > !!$             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
877 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
878 > !!$        
879 > !!$          else
880               pref =  pre14 * q_i / 3.0_dp
881               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
882                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
# Line 865 | Line 909 | contains
909               duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
910               duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
911            
912 <          endif
912 > !!$          endif
913         endif
914      endif
915  
# Line 875 | Line 919 | contains
919            
920            pref = pre12 * q_j * mu_i
921            
922 <          if (summationMethod .eq. UNDAMPED_WOLF) then
922 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
923 > !!$             ri2 = riji * riji
924 > !!$             ri3 = ri2 * riji
925 > !!$
926 > !!$             pref = pre12 * q_j * mu_i
927 > !!$             vterm = pref * ct_i * (ri2 - rcuti2)
928 > !!$             vpair = vpair + vterm
929 > !!$             epot = epot + sw*vterm
930 > !!$            
931 > !!$             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
932 > !!$                  - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
933 > !!$             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
934 > !!$                  - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
935 > !!$             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
936 > !!$                  - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
937 > !!$            
938 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*( ri2*xhat - d(1)*rcuti3 )
939 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*( ri2*yhat - d(2)*rcuti3 )
940 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*( ri2*zhat - d(3)*rcuti3 )
941 > !!$
942 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
943 >          if (summationMethod .eq. REACTION_FIELD) then
944               ri2 = riji * riji
945               ri3 = ri2 * riji
946  
947               pref = pre12 * q_j * mu_i
883 –             vterm = pref * ct_i * (ri2 - rcuti2)
884 –             vpair = vpair + vterm
885 –             epot = epot + sw*vterm
886 –            
887 –             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
888 –                  - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
889 –             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
890 –                  - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
891 –             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
892 –                  - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
893 –            
894 –             duduz_i(1) = duduz_i(1) + sw*pref*( ri2*xhat - d(1)*rcuti3 )
895 –             duduz_i(2) = duduz_i(2) + sw*pref*( ri2*yhat - d(2)*rcuti3 )
896 –             duduz_i(3) = duduz_i(3) + sw*pref*( ri2*zhat - d(3)*rcuti3 )
897 –
898 –          elseif (summationMethod .eq. REACTION_FIELD) then
899 –             ri2 = riji * riji
900 –             ri3 = ri2 * riji
901 –
902 –             pref = pre12 * q_j * mu_i
948               vterm = pref * ct_i * ( ri2 - preRF2*rij )
949               vpair = vpair + vterm
950               epot = epot + sw*vterm
# Line 946 | Line 991 | contains
991        
992         if (j_is_Dipole) then
993  
994 <          if (summationMethod .eq. UNDAMPED_WOLF) then
994 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
995   !!$             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
996   !!$
997   !!$             ri2 = riji * riji
# Line 982 | Line 1027 | contains
1027   !!$                  - rcuti3*(uz_i(2) - 3.0d0*ct_i*yhat))
1028   !!$             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1029   !!$                  - rcuti3*(uz_i(3) - 3.0d0*ct_i*zhat))
1030 <          
1031 <          elseif (summationMethod .eq. DAMPED_WOLF) then
1032 <             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1033 <            
1034 <             ri2 = riji * riji
1035 <             ri3 = ri2 * riji
1036 <             ri4 = ri2 * ri2
1037 <             sc2 = scale * scale
1038 <            
1039 <             pref = pre22 * mu_i * mu_j
1040 <             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j)
1041 <             vpair = vpair + vterm
1042 <             epot = epot + sw*vterm
1043 <            
1044 <             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1045 <            
1046 <             dudx = dudx + sw*pref*3.0d0*ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1047 <             dudy = dudy + sw*pref*3.0d0*ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1048 <             dudz = dudz + sw*pref*3.0d0*ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1049 <            
1050 <             duduz_i(1) = duduz_i(1) + sw*pref*ri3 *(uz_j(1) - 3.0d0*ct_j*xhat)
1051 <             duduz_i(2) = duduz_i(2) + sw*pref*ri3 *(uz_j(2) - 3.0d0*ct_j*yhat)
1052 <             duduz_i(3) = duduz_i(3) + sw*pref*ri3 *(uz_j(3) - 3.0d0*ct_j*zhat)
1053 <            
1054 <             duduz_j(1) = duduz_j(1) + sw*pref*ri3 *(uz_i(1) - 3.0d0*ct_i*xhat)
1055 <             duduz_j(2) = duduz_j(2) + sw*pref*ri3 *(uz_i(2) - 3.0d0*ct_i*yhat)
1056 <             duduz_j(3) = duduz_j(3) + sw*pref*ri3 *(uz_i(3) - 3.0d0*ct_i*zhat)
1057 <            
1058 <          elseif (summationMethod .eq. REACTION_FIELD) then
1030 > !!$          
1031 > !!$          elseif (summationMethod .eq. DAMPED_WOLF) then
1032 > !!$             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1033 > !!$            
1034 > !!$             ri2 = riji * riji
1035 > !!$             ri3 = ri2 * riji
1036 > !!$             ri4 = ri2 * ri2
1037 > !!$             sc2 = scale * scale
1038 > !!$            
1039 > !!$             pref = pre22 * mu_i * mu_j
1040 > !!$             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j)
1041 > !!$             vpair = vpair + vterm
1042 > !!$             epot = epot + sw*vterm
1043 > !!$            
1044 > !!$             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1045 > !!$            
1046 > !!$             dudx = dudx + sw*pref*3.0d0*ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1047 > !!$             dudy = dudy + sw*pref*3.0d0*ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1048 > !!$             dudz = dudz + sw*pref*3.0d0*ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1049 > !!$            
1050 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*ri3 *(uz_j(1) - 3.0d0*ct_j*xhat)
1051 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*ri3 *(uz_j(2) - 3.0d0*ct_j*yhat)
1052 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*ri3 *(uz_j(3) - 3.0d0*ct_j*zhat)
1053 > !!$            
1054 > !!$             duduz_j(1) = duduz_j(1) + sw*pref*ri3 *(uz_i(1) - 3.0d0*ct_i*xhat)
1055 > !!$             duduz_j(2) = duduz_j(2) + sw*pref*ri3 *(uz_i(2) - 3.0d0*ct_i*yhat)
1056 > !!$             duduz_j(3) = duduz_j(3) + sw*pref*ri3 *(uz_i(3) - 3.0d0*ct_i*zhat)
1057 > !!$            
1058 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
1059 >          if (summationMethod .eq. REACTION_FIELD) then
1060               ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1061  
1062               ri2 = riji * riji
# Line 1114 | Line 1160 | contains
1160            cy2 = cy_i * cy_i
1161            cz2 = cz_i * cz_i
1162  
1163 <          if (summationMethod .eq. UNDAMPED_WOLF) then
1164 <             pref = pre14 * q_j / 3.0_dp
1165 <             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1166 <                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1167 <                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1168 <             vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1169 <                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1170 <                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1171 <             vpair = vpair + ( vterm1 - vterm2 )
1172 <             epot = epot + sw*( vterm1 - vterm2 )
1173 <            
1174 <             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1175 <                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1176 <                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1177 <                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1178 <                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1179 <                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1180 <                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1181 <             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1182 <                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1183 <                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1184 <                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1185 <                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1186 <                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1187 <                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1188 <             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1189 <                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1190 <                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1191 <                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1192 <                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1193 <                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1194 <                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1195 <            
1196 <             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1197 <                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1198 <             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1199 <                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1200 <             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1201 <                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1202 <            
1203 <             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1204 <                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1205 <             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1206 <                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1207 <             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1208 <                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1209 <            
1210 <             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1211 <                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1212 <             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1213 <                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1214 <             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1215 <                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1216 <
1217 <          else
1163 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
1164 > !!$             pref = pre14 * q_j / 3.0_dp
1165 > !!$             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1166 > !!$                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1167 > !!$                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1168 > !!$             vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1169 > !!$                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1170 > !!$                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1171 > !!$             vpair = vpair + ( vterm1 - vterm2 )
1172 > !!$             epot = epot + sw*( vterm1 - vterm2 )
1173 > !!$            
1174 > !!$             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1175 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1176 > !!$                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1177 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1178 > !!$                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1179 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1180 > !!$                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1181 > !!$             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1182 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1183 > !!$                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1184 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1185 > !!$                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1186 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1187 > !!$                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1188 > !!$             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1189 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1190 > !!$                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1191 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1192 > !!$                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1193 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1194 > !!$                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1195 > !!$            
1196 > !!$             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1197 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1198 > !!$             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1199 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1200 > !!$             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1201 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1202 > !!$            
1203 > !!$             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1204 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1205 > !!$             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1206 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1207 > !!$             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1208 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1209 > !!$            
1210 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1211 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1212 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1213 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1214 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1215 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1216 > !!$
1217 > !!$          else
1218               pref = pre14 * q_j / 3.0_dp
1219               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1220                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
# Line 1200 | Line 1246 | contains
1246               duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1247               duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1248               duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1249 <          endif
1249 > !!$          endif
1250         endif
1251      endif
1252  
# Line 1361 | Line 1407 | contains
1407            
1408         endif
1409  
1410 <    elseif (summationMethod .eq. UNDAMPED_WOLF) then
1411 <       if (ElectrostaticMap(atid1)%is_Charge) then
1412 <          c1 = getCharge(atid1)
1413 <          
1414 <          mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1415 <       endif
1416 <      
1417 <    elseif (summationMethod .eq. DAMPED_WOLF) then
1418 <       if (ElectrostaticMap(atid1)%is_Charge) then
1419 <          c1 = getCharge(atid1)
1420 <          
1421 <          mypot = mypot - (constERFC * rcuti * 0.5_dp + &
1422 <               dampingAlpha*invRootPi) * c1 * c1      
1423 <       endif
1410 > !!$    elseif (summationMethod .eq. UNDAMPED_WOLF) then
1411 > !!$       if (ElectrostaticMap(atid1)%is_Charge) then
1412 > !!$          c1 = getCharge(atid1)
1413 > !!$          
1414 > !!$          mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1415 > !!$       endif
1416 > !!$      
1417 > !!$    elseif (summationMethod .eq. DAMPED_WOLF) then
1418 > !!$       if (ElectrostaticMap(atid1)%is_Charge) then
1419 > !!$          c1 = getCharge(atid1)
1420 > !!$          
1421 > !!$          mypot = mypot - (constERFC * rcuti * 0.5_dp + &
1422 > !!$               dampingAlpha*invRootPi) * c1 * c1      
1423 > !!$       endif
1424      endif
1425      
1426      return

Diff Legend

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