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 2405 by chrisfen, Tue Nov 1 19:24:57 2005 UTC vs.
Revision 2409 by chrisfen, Wed Nov 2 20:36:25 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 76 | Line 77 | module electrostatic_module
77  
78    !! variables to handle different summation methods for long-range electrostatics:
79    integer, save :: summationMethod = NONE
80 +  integer, save :: screeningMethod = UNDAMPED
81    logical, save :: summationMethodChecked = .false.
82    real(kind=DP), save :: defaultCutoff = 0.0_DP
83    real(kind=DP), save :: defaultCutoff2 = 0.0_DP
# Line 104 | Line 106 | module electrostatic_module
106   #endif
107    
108    public :: setElectrostaticSummationMethod
109 +  public :: setScreeningMethod
110    public :: setElectrostaticCutoffRadius
111 <  public :: setDampedWolfAlpha
111 >  public :: setDampingAlpha
112    public :: setReactionFieldDielectric
113    public :: newElectrostaticType
114    public :: setCharge
# Line 146 | Line 149 | contains
149      summationMethod = the_ESM
150  
151    end subroutine setElectrostaticSummationMethod
152 +
153 +  subroutine setScreeningMethod(the_SM)
154 +    integer, intent(in) :: the_SM    
155 +    screeningMethod = the_SM
156 +  end subroutine setScreeningMethod
157  
158    subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
159      real(kind=dp), intent(in) :: thisRcut
# Line 156 | Line 164 | contains
164      haveDefaultCutoff = .true.
165    end subroutine setElectrostaticCutoffRadius
166  
167 <  subroutine setDampedWolfAlpha(thisAlpha)
167 >  subroutine setDampingAlpha(thisAlpha)
168      real(kind=dp), intent(in) :: thisAlpha
169      dampingAlpha = thisAlpha
170      haveDampingAlpha = .true.
171 <  end subroutine setDampedWolfAlpha
171 >  end subroutine setDampingAlpha
172    
173    subroutine setReactionFieldDielectric(thisDielectric)
174      real(kind=dp), intent(in) :: thisDielectric
# Line 389 | Line 397 | contains
397      rcuti3 = rcuti2*rcuti
398      rcuti4 = rcuti2*rcuti2
399  
400 <    if (summationMethod .eq. DAMPED_WOLF) then
400 >    if (screeningMethod .eq. DAMPED) then
401         if (.not.haveDampingAlpha) then
402            call handleError("checkSummationMethod", "no Damping Alpha set!")
403         endif
# Line 420 | Line 428 | contains
428      summationMethodChecked = .true.
429    end subroutine checkSummationMethod
430  
431 <
431 > !!$
432 > !!$  subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
433 > !!$       vpair, fpair, pot, eFrame, f, t, do_pot)
434    subroutine doElectrostaticPair(atom1, atom2, d, rij, r2, sw, &
435 <       vpair, fpair, pot, eFrame, f, t, do_pot)
435 >       vpair, fpair, pot, eFrame, f, t, do_pot, fstrs)
436  
437      logical, intent(in) :: do_pot
438  
# Line 433 | Line 443 | contains
443      real(kind=dp), intent(in), dimension(3) :: d
444      real(kind=dp), intent(inout) :: vpair
445      real(kind=dp), intent(inout), dimension(3) :: fpair    
446 +    real(kind=dp), intent(inout), dimension(3) :: fstrs
447  
448      real( kind = dp ) :: pot
449      real( kind = dp ), dimension(9,nLocal) :: eFrame
450      real( kind = dp ), dimension(3,nLocal) :: f
451 +    real( kind = dp ), dimension(3,nLocal) :: felec
452      real( kind = dp ), dimension(3,nLocal) :: t
453  
454      real (kind = dp), dimension(3) :: ux_i, uy_i, uz_i
# Line 481 | Line 493 | contains
493      me2 = atid(atom2)
494   #endif
495  
496 + !!$    if (rij .ge. defaultCutoff) then
497 + !!$       write(*,*) 'warning: rij = ', rij, ' rcut = ', defaultCutoff, ' sw = ', sw
498 + !!$    endif
499 +
500      !! some variables we'll need independent of electrostatic type:
501  
502      riji = 1.0d0 / rij
# Line 624 | Line 640 | contains
640  
641         if (j_is_Charge) then
642  
643 <          if (summationMethod .eq. UNDAMPED_WOLF) then
643 >          if (summationMethod .eq. SHIFTED_POTENTIAL) then
644 >
645               vterm = pre11 * q_i * q_j * (riji - rcuti)
646               vpair = vpair + vterm
647               epot = epot + sw*vterm
648              
649 <             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)
649 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji)
650              
651               dudx = dudx + dudr * xhat
652               dudy = dudy + dudr * yhat
653               dudz = dudz + dudr * zhat
654  
655 <          elseif (summationMethod .eq. DAMPED_WOLF) then
656 <             varERFC = derfc(dampingAlpha*rij)
657 <             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
658 <             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
655 > !!$          elseif (summationMethod .eq. DAMPED_WOLF) then
656 > !!$             varERFC = derfc(dampingAlpha*rij)
657 > !!$             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
658 > !!$             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
659 > !!$             vpair = vpair + vterm
660 > !!$             epot = epot + sw*vterm
661 > !!$            
662 > !!$!             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
663 > !!$!                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
664 > !!$!                  + alphaPi*constEXP*rcuti)) )
665 > !!$             dudr  = -sw*pre11*q_i*q_j * (varERFC*riji*riji &
666 > !!$                  + alphaPi*varEXP*riji)
667 > !!$            
668 > !!$             dudx = dudx + dudr * xhat
669 > !!$             dudy = dudy + dudr * yhat
670 > !!$             dudz = dudz + dudr * zhat
671 >
672 >          elseif (summationMethod .eq. SHIFTED_FORCE) then
673 >             vterm = pre11 * q_i * q_j * (riji + rij*rcuti2 - 2.0d0*rcuti)
674               vpair = vpair + vterm
675               epot = epot + sw*vterm
676              
677 <             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
678 <                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
647 <                  + alphaPi*constEXP*rcuti)) )
648 <            
677 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)
678 >                          
679               dudx = dudx + dudr * xhat
680               dudy = dudy + dudr * yhat
681               dudz = dudz + dudr * zhat
682  
683 + !!$          elseif (summationMethod .eq. DAMPED_SHIFTED_FORCE) then
684 + !!$             varERFC = derfc(dampingAlpha*rij)
685 + !!$             varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
686 + !!$             vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
687 + !!$             vpair = vpair + vterm
688 + !!$             epot = epot + sw*vterm
689 + !!$            
690 + !!$!             dudr  = -sw*pre11*q_i*q_j * (((varERFC*riji*riji &
691 + !!$!                  + alphaPi*varEXP*riji) - (constERFC*rcuti2 &
692 + !!$!                  + alphaPi*constEXP*rcuti)) )
693 + !!$             dudr  = -sw*pre11*q_i*q_j * (varERFC*riji*riji &
694 + !!$                  + alphaPi*varEXP*riji)
695 + !!$            
696 + !!$             dudx = dudx + dudr * xhat
697 + !!$             dudy = dudy + dudr * yhat
698 + !!$             dudz = dudz + dudr * zhat
699 +
700            elseif (summationMethod .eq. REACTION_FIELD) then
701               preVal = pre11 * q_i * q_j
702               rfVal = preRF*rij*rij
# Line 683 | Line 730 | contains
730  
731            pref = pre12 * q_i * mu_j
732  
733 <          if (summationMethod .eq. UNDAMPED_WOLF) then
734 <             ri2 = riji * riji
735 <             ri3 = ri2 * riji
736 <
737 <             pref = pre12 * q_i * mu_j
738 <             vterm = - pref * ct_j * (ri2 - rcuti2)
739 <             vpair = vpair + vterm
740 <             epot = epot + sw*vterm
741 <            
742 <             !! this has a + sign in the () because the rij vector is
743 <             !! r_j - r_i and the charge-dipole potential takes the origin
744 <             !! as the point dipole, which is atom j in this case.
745 <            
746 <             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
747 <                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
748 <             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
749 <                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
750 <             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
751 <                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
752 <            
753 <             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
754 <             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
755 <             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
733 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
734 > !!$             ri2 = riji * riji
735 > !!$             ri3 = ri2 * riji
736 > !!$
737 > !!$             pref = pre12 * q_i * mu_j
738 > !!$             vterm = - pref * ct_j * (ri2 - rcuti2)
739 > !!$             vpair = vpair + vterm
740 > !!$             epot = epot + sw*vterm
741 > !!$            
742 > !!$             !! this has a + sign in the () because the rij vector is
743 > !!$             !! r_j - r_i and the charge-dipole potential takes the origin
744 > !!$             !! as the point dipole, which is atom j in this case.
745 > !!$            
746 > !!$             dudx = dudx - sw*pref * ( ri3*( uz_j(1) - 3.0d0*ct_j*xhat) &
747 > !!$                  - rcuti3*( uz_j(1) - 3.0d0*ct_j*d(1)*rcuti ) )
748 > !!$             dudy = dudy - sw*pref * ( ri3*( uz_j(2) - 3.0d0*ct_j*yhat) &
749 > !!$                  - rcuti3*( uz_j(2) - 3.0d0*ct_j*d(2)*rcuti ) )
750 > !!$             dudz = dudz - sw*pref * ( ri3*( uz_j(3) - 3.0d0*ct_j*zhat) &
751 > !!$                  - rcuti3*( uz_j(3) - 3.0d0*ct_j*d(3)*rcuti ) )
752 > !!$            
753 > !!$             duduz_j(1) = duduz_j(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
754 > !!$             duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
755 > !!$             duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
756 > !!$
757 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
758  
759 <          elseif (summationMethod .eq. REACTION_FIELD) then
759 >          if (summationMethod .eq. REACTION_FIELD) then
760               ri2 = riji * riji
761               ri3 = ri2 * riji
762      
# Line 772 | Line 821 | contains
821            cy2 = cy_j * cy_j
822            cz2 = cz_j * cz_j
823  
824 <          if (summationMethod .eq. UNDAMPED_WOLF) then
824 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
825 > !!$             pref =  pre14 * q_i / 3.0_dp
826 > !!$             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
827 > !!$                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
828 > !!$                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
829 > !!$             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
830 > !!$                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
831 > !!$                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
832 > !!$             vpair = vpair + ( vterm1 - vterm2 )
833 > !!$             epot = epot + sw*( vterm1 - vterm2 )
834 > !!$            
835 > !!$             dudx = dudx - (5.0_dp * &
836 > !!$                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
837 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
838 > !!$                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
839 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
840 > !!$                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
841 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
842 > !!$                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
843 > !!$             dudy = dudy - (5.0_dp * &
844 > !!$                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
845 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
846 > !!$                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
847 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
848 > !!$                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
849 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
850 > !!$                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
851 > !!$             dudz = dudz - (5.0_dp * &
852 > !!$                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
853 > !!$                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
854 > !!$                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
855 > !!$                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
856 > !!$                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
857 > !!$                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
858 > !!$                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
859 > !!$            
860 > !!$             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
861 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
862 > !!$             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
863 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
864 > !!$             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
865 > !!$                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
866 > !!$            
867 > !!$             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
868 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
869 > !!$             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
870 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
871 > !!$             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
872 > !!$                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
873 > !!$            
874 > !!$             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
875 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
876 > !!$             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
877 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
878 > !!$             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
879 > !!$                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
880 > !!$        
881 > !!$          else
882               pref =  pre14 * q_i / 3.0_dp
777 –             vterm1 = pref * ri3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
778 –                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
779 –                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
780 –             vterm2 = pref * rcuti3*( qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
781 –                  qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
782 –                  qzz_j * (3.0_dp*cz2 - 1.0_dp) )
783 –             vpair = vpair + ( vterm1 - vterm2 )
784 –             epot = epot + sw*( vterm1 - vterm2 )
785 –            
786 –             dudx = dudx - (5.0_dp * &
787 –                  (vterm1*riji*xhat - vterm2*rcuti2*d(1))) + sw*pref * ( &
788 –                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(1)) - &
789 –                  qxx_j*2.0_dp*(xhat - rcuti*d(1))) + &
790 –                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(1)) - &
791 –                  qyy_j*2.0_dp*(xhat - rcuti*d(1))) + &
792 –                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(1)) - &
793 –                  qzz_j*2.0_dp*(xhat - rcuti*d(1))) )
794 –             dudy = dudy - (5.0_dp * &
795 –                  (vterm1*riji*yhat - vterm2*rcuti2*d(2))) + sw*pref * ( &
796 –                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(2)) - &
797 –                  qxx_j*2.0_dp*(yhat - rcuti*d(2))) + &
798 –                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(2)) - &
799 –                  qyy_j*2.0_dp*(yhat - rcuti*d(2))) + &
800 –                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(2)) - &
801 –                  qzz_j*2.0_dp*(yhat - rcuti*d(2))) )
802 –             dudz = dudz - (5.0_dp * &
803 –                  (vterm1*riji*zhat - vterm2*rcuti2*d(3))) + sw*pref * ( &
804 –                  (ri4 - rcuti4)*(qxx_j*(6.0_dp*cx_j*ux_j(3)) - &
805 –                  qxx_j*2.0_dp*(zhat - rcuti*d(3))) + &
806 –                  (ri4 - rcuti4)*(qyy_j*(6.0_dp*cy_j*uy_j(3)) - &
807 –                  qyy_j*2.0_dp*(zhat - rcuti*d(3))) + &
808 –                  (ri4 - rcuti4)*(qzz_j*(6.0_dp*cz_j*uz_j(3)) - &
809 –                  qzz_j*2.0_dp*(zhat - rcuti*d(3))) )
810 –            
811 –             dudux_j(1) = dudux_j(1) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*xhat) -&
812 –                  rcuti4*(qxx_j*6.0_dp*cx_j*d(1)))
813 –             dudux_j(2) = dudux_j(2) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*yhat) -&
814 –                  rcuti4*(qxx_j*6.0_dp*cx_j*d(2)))
815 –             dudux_j(3) = dudux_j(3) + sw*pref*(ri3*(qxx_j*6.0_dp*cx_j*zhat) -&
816 –                  rcuti4*(qxx_j*6.0_dp*cx_j*d(3)))
817 –            
818 –             duduy_j(1) = duduy_j(1) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*xhat) -&
819 –                  rcuti4*(qyy_j*6.0_dp*cx_j*d(1)))
820 –             duduy_j(2) = duduy_j(2) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*yhat) -&
821 –                  rcuti4*(qyy_j*6.0_dp*cx_j*d(2)))
822 –             duduy_j(3) = duduy_j(3) + sw*pref*(ri3*(qyy_j*6.0_dp*cy_j*zhat) -&
823 –                  rcuti4*(qyy_j*6.0_dp*cx_j*d(3)))
824 –            
825 –             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*xhat) -&
826 –                  rcuti4*(qzz_j*6.0_dp*cx_j*d(1)))
827 –             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*yhat) -&
828 –                  rcuti4*(qzz_j*6.0_dp*cx_j*d(2)))
829 –             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(qzz_j*6.0_dp*cz_j*zhat) -&
830 –                  rcuti4*(qzz_j*6.0_dp*cx_j*d(3)))
831 –        
832 –          else
833 –             pref =  pre14 * q_i / 3.0_dp
883               vterm = pref * ri3 * (qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
884                    qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
885                    qzz_j * (3.0_dp*cz2 - 1.0_dp))
# Line 862 | Line 911 | contains
911               duduz_j(2) = duduz_j(2) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*yhat)
912               duduz_j(3) = duduz_j(3) + sw*pref * ri3*(qzz_j*6.0_dp*cz_j*zhat)
913            
914 <          endif
914 > !!$          endif
915         endif
916      endif
917  
# Line 872 | Line 921 | contains
921            
922            pref = pre12 * q_j * mu_i
923            
924 <          if (summationMethod .eq. UNDAMPED_WOLF) then
924 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
925 > !!$             ri2 = riji * riji
926 > !!$             ri3 = ri2 * riji
927 > !!$
928 > !!$             pref = pre12 * q_j * mu_i
929 > !!$             vterm = pref * ct_i * (ri2 - rcuti2)
930 > !!$             vpair = vpair + vterm
931 > !!$             epot = epot + sw*vterm
932 > !!$            
933 > !!$             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
934 > !!$                  - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
935 > !!$             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
936 > !!$                  - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
937 > !!$             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
938 > !!$                  - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
939 > !!$            
940 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*( ri2*xhat - d(1)*rcuti3 )
941 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*( ri2*yhat - d(2)*rcuti3 )
942 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*( ri2*zhat - d(3)*rcuti3 )
943 > !!$
944 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
945 >          if (summationMethod .eq. REACTION_FIELD) then
946               ri2 = riji * riji
947               ri3 = ri2 * riji
948  
949               pref = pre12 * q_j * mu_i
880 –             vterm = pref * ct_i * (ri2 - rcuti2)
881 –             vpair = vpair + vterm
882 –             epot = epot + sw*vterm
883 –            
884 –             dudx = dudx + sw*pref * ( ri3*( uz_i(1) - 3.0d0*ct_i*xhat) &
885 –                  - rcuti3*( uz_i(1) - 3.0d0*ct_i*d(1)*rcuti ) )
886 –             dudy = dudy + sw*pref * ( ri3*( uz_i(2) - 3.0d0*ct_i*yhat) &
887 –                  - rcuti3*( uz_i(2) - 3.0d0*ct_i*d(2)*rcuti ) )
888 –             dudz = dudz + sw*pref * ( ri3*( uz_i(3) - 3.0d0*ct_i*zhat) &
889 –                  - rcuti3*( uz_i(3) - 3.0d0*ct_i*d(3)*rcuti ) )
890 –            
891 –             duduz_i(1) = duduz_i(1) + sw*pref*( ri2*xhat - d(1)*rcuti3 )
892 –             duduz_i(2) = duduz_i(2) + sw*pref*( ri2*yhat - d(2)*rcuti3 )
893 –             duduz_i(3) = duduz_i(3) + sw*pref*( ri2*zhat - d(3)*rcuti3 )
894 –
895 –          elseif (summationMethod .eq. REACTION_FIELD) then
896 –             ri2 = riji * riji
897 –             ri3 = ri2 * riji
898 –
899 –             pref = pre12 * q_j * mu_i
950               vterm = pref * ct_i * ( ri2 - preRF2*rij )
951               vpair = vpair + vterm
952               epot = epot + sw*vterm
# Line 943 | Line 993 | contains
993        
994         if (j_is_Dipole) then
995  
996 <          if (summationMethod .eq. UNDAMPED_WOLF) then
996 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
997   !!$             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
998   !!$
999   !!$             ri2 = riji * riji
# Line 979 | Line 1029 | contains
1029   !!$                  - rcuti3*(uz_i(2) - 3.0d0*ct_i*yhat))
1030   !!$             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1031   !!$                  - rcuti3*(uz_i(3) - 3.0d0*ct_i*zhat))
1032 <          
1033 <          elseif (summationMethod .eq. DAMPED_WOLF) then
1034 <             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1035 <            
1036 <             ri2 = riji * riji
1037 <             ri3 = ri2 * riji
1038 <             ri4 = ri2 * ri2
1039 <             sc2 = scale * scale
1040 <            
1041 <             pref = pre22 * mu_i * mu_j
1042 <             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j)
1043 <             vpair = vpair + vterm
1044 <             epot = epot + sw*vterm
1045 <            
1046 <             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1047 <            
1048 <             dudx = dudx + sw*pref*3.0d0*ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1049 <             dudy = dudy + sw*pref*3.0d0*ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1050 <             dudz = dudz + sw*pref*3.0d0*ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1051 <            
1052 <             duduz_i(1) = duduz_i(1) + sw*pref*ri3 *(uz_j(1) - 3.0d0*ct_j*xhat)
1053 <             duduz_i(2) = duduz_i(2) + sw*pref*ri3 *(uz_j(2) - 3.0d0*ct_j*yhat)
1054 <             duduz_i(3) = duduz_i(3) + sw*pref*ri3 *(uz_j(3) - 3.0d0*ct_j*zhat)
1055 <            
1056 <             duduz_j(1) = duduz_j(1) + sw*pref*ri3 *(uz_i(1) - 3.0d0*ct_i*xhat)
1057 <             duduz_j(2) = duduz_j(2) + sw*pref*ri3 *(uz_i(2) - 3.0d0*ct_i*yhat)
1058 <             duduz_j(3) = duduz_j(3) + sw*pref*ri3 *(uz_i(3) - 3.0d0*ct_i*zhat)
1059 <            
1060 <          elseif (summationMethod .eq. REACTION_FIELD) then
1032 > !!$          
1033 > !!$          elseif (summationMethod .eq. DAMPED_WOLF) then
1034 > !!$             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1035 > !!$            
1036 > !!$             ri2 = riji * riji
1037 > !!$             ri3 = ri2 * riji
1038 > !!$             ri4 = ri2 * ri2
1039 > !!$             sc2 = scale * scale
1040 > !!$            
1041 > !!$             pref = pre22 * mu_i * mu_j
1042 > !!$             vterm = pref * ri3 * (ct_ij - 3.0d0 * ct_i * ct_j)
1043 > !!$             vpair = vpair + vterm
1044 > !!$             epot = epot + sw*vterm
1045 > !!$            
1046 > !!$             a1 = 5.0d0 * ct_i * ct_j * sc2 - ct_ij
1047 > !!$            
1048 > !!$             dudx = dudx + sw*pref*3.0d0*ri4*(a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1049 > !!$             dudy = dudy + sw*pref*3.0d0*ri4*(a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1050 > !!$             dudz = dudz + sw*pref*3.0d0*ri4*(a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1051 > !!$            
1052 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*ri3 *(uz_j(1) - 3.0d0*ct_j*xhat)
1053 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*ri3 *(uz_j(2) - 3.0d0*ct_j*yhat)
1054 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*ri3 *(uz_j(3) - 3.0d0*ct_j*zhat)
1055 > !!$            
1056 > !!$             duduz_j(1) = duduz_j(1) + sw*pref*ri3 *(uz_i(1) - 3.0d0*ct_i*xhat)
1057 > !!$             duduz_j(2) = duduz_j(2) + sw*pref*ri3 *(uz_i(2) - 3.0d0*ct_i*yhat)
1058 > !!$             duduz_j(3) = duduz_j(3) + sw*pref*ri3 *(uz_i(3) - 3.0d0*ct_i*zhat)
1059 > !!$            
1060 > !!$          elseif (summationMethod .eq. REACTION_FIELD) then
1061 >          if (summationMethod .eq. REACTION_FIELD) then
1062               ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1063  
1064               ri2 = riji * riji
# Line 1111 | Line 1162 | contains
1162            cy2 = cy_i * cy_i
1163            cz2 = cz_i * cz_i
1164  
1165 <          if (summationMethod .eq. UNDAMPED_WOLF) then
1166 <             pref = pre14 * q_j / 3.0_dp
1167 <             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1168 <                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1169 <                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1170 <             vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1171 <                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1172 <                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1173 <             vpair = vpair + ( vterm1 - vterm2 )
1174 <             epot = epot + sw*( vterm1 - vterm2 )
1175 <            
1176 <             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1177 <                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1178 <                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1179 <                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1180 <                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1181 <                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1182 <                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1183 <             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1184 <                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1185 <                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1186 <                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1187 <                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1188 <                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1189 <                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1190 <             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1191 <                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1192 <                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1193 <                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1194 <                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1195 <                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1196 <                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1197 <            
1198 <             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1199 <                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1200 <             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1201 <                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1202 <             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1203 <                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1204 <            
1205 <             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1206 <                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1207 <             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1208 <                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1209 <             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1210 <                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1211 <            
1212 <             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1213 <                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1214 <             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1215 <                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1216 <             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1217 <                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1218 <
1219 <          else
1165 > !!$          if (summationMethod .eq. UNDAMPED_WOLF) then
1166 > !!$             pref = pre14 * q_j / 3.0_dp
1167 > !!$             vterm1 = pref * ri3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1168 > !!$                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1169 > !!$                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1170 > !!$             vterm2 = pref * rcuti3*( qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1171 > !!$                  qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1172 > !!$                  qzz_i * (3.0_dp*cz2 - 1.0_dp) )
1173 > !!$             vpair = vpair + ( vterm1 - vterm2 )
1174 > !!$             epot = epot + sw*( vterm1 - vterm2 )
1175 > !!$            
1176 > !!$             dudx = dudx - sw*(5.0_dp*(vterm1*riji*xhat-vterm2*rcuti2*d(1))) +&
1177 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(1)) - &
1178 > !!$                  qxx_i*2.0_dp*(xhat - rcuti*d(1))) + &
1179 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(1)) - &
1180 > !!$                  qyy_i*2.0_dp*(xhat - rcuti*d(1))) + &
1181 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(1)) - &
1182 > !!$                  qzz_i*2.0_dp*(xhat - rcuti*d(1))) )
1183 > !!$             dudy = dudy - sw*(5.0_dp*(vterm1*riji*yhat-vterm2*rcuti2*d(2))) +&
1184 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(2)) - &
1185 > !!$                  qxx_i*2.0_dp*(yhat - rcuti*d(2))) + &
1186 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(2)) - &
1187 > !!$                  qyy_i*2.0_dp*(yhat - rcuti*d(2))) + &
1188 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(2)) - &
1189 > !!$                  qzz_i*2.0_dp*(yhat - rcuti*d(2))) )
1190 > !!$             dudz = dudz - sw*(5.0_dp*(vterm1*riji*zhat-vterm2*rcuti2*d(3))) +&
1191 > !!$                  sw*pref * ( (ri4 - rcuti4)*(qxx_i*(6.0_dp*cx_i*ux_i(3)) - &
1192 > !!$                  qxx_i*2.0_dp*(zhat - rcuti*d(3))) + &
1193 > !!$                  (ri4 - rcuti4)*(qyy_i*(6.0_dp*cy_i*uy_i(3)) - &
1194 > !!$                  qyy_i*2.0_dp*(zhat - rcuti*d(3))) + &
1195 > !!$                  (ri4 - rcuti4)*(qzz_i*(6.0_dp*cz_i*uz_i(3)) - &
1196 > !!$                  qzz_i*2.0_dp*(zhat - rcuti*d(3))) )
1197 > !!$            
1198 > !!$             dudux_i(1) = dudux_i(1) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*xhat) -&
1199 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(1)))
1200 > !!$             dudux_i(2) = dudux_i(2) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*yhat) -&
1201 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(2)))
1202 > !!$             dudux_i(3) = dudux_i(3) + sw*pref*(ri3*(qxx_i*6.0_dp*cx_i*zhat) -&
1203 > !!$                  rcuti4*(qxx_i*6.0_dp*cx_i*d(3)))
1204 > !!$            
1205 > !!$             duduy_i(1) = duduy_i(1) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*xhat) -&
1206 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(1)))
1207 > !!$             duduy_i(2) = duduy_i(2) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*yhat) -&
1208 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(2)))
1209 > !!$             duduy_i(3) = duduy_i(3) + sw*pref*(ri3*(qyy_i*6.0_dp*cy_i*zhat) -&
1210 > !!$                  rcuti4*(qyy_i*6.0_dp*cx_i*d(3)))
1211 > !!$            
1212 > !!$             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*xhat) -&
1213 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(1)))
1214 > !!$             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*yhat) -&
1215 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(2)))
1216 > !!$             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(qzz_i*6.0_dp*cz_i*zhat) -&
1217 > !!$                  rcuti4*(qzz_i*6.0_dp*cx_i*d(3)))
1218 > !!$
1219 > !!$          else
1220               pref = pre14 * q_j / 3.0_dp
1221               vterm = pref * ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1222                    qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
# Line 1197 | Line 1248 | contains
1248               duduz_i(1) = duduz_i(1) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*xhat)
1249               duduz_i(2) = duduz_i(2) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*yhat)
1250               duduz_i(3) = duduz_i(3) + sw*pref*ri3*(qzz_i*6.0_dp*cz_i*zhat)
1251 <          endif
1251 > !!$          endif
1252         endif
1253      endif
1254  
# Line 1358 | Line 1409 | contains
1409            
1410         endif
1411  
1412 <    elseif (summationMethod .eq. UNDAMPED_WOLF) then
1413 <       if (ElectrostaticMap(atid1)%is_Charge) then
1414 <          c1 = getCharge(atid1)
1415 <          
1416 <          mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1417 <       endif
1418 <      
1419 <    elseif (summationMethod .eq. DAMPED_WOLF) then
1420 <       if (ElectrostaticMap(atid1)%is_Charge) then
1421 <          c1 = getCharge(atid1)
1422 <          
1423 <          mypot = mypot - (constERFC * rcuti * 0.5_dp + &
1424 <               dampingAlpha*invRootPi) * c1 * c1      
1425 <       endif
1412 > !!$    elseif (summationMethod .eq. UNDAMPED_WOLF) then
1413 > !!$       if (ElectrostaticMap(atid1)%is_Charge) then
1414 > !!$          c1 = getCharge(atid1)
1415 > !!$          
1416 > !!$          mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1417 > !!$       endif
1418 > !!$      
1419 > !!$    elseif (summationMethod .eq. DAMPED_WOLF) then
1420 > !!$       if (ElectrostaticMap(atid1)%is_Charge) then
1421 > !!$          c1 = getCharge(atid1)
1422 > !!$          
1423 > !!$          mypot = mypot - (constERFC * rcuti * 0.5_dp + &
1424 > !!$               dampingAlpha*invRootPi) * c1 * c1      
1425 > !!$       endif
1426      endif
1427      
1428      return

Diff Legend

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