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 2552 by chrisfen, Thu Jan 12 16:47:25 2006 UTC vs.
Revision 2717 by gezelter, Mon Apr 17 21:49:12 2006 UTC

# Line 47 | Line 47 | module electrostatic_module
47    use vector_class
48    use simulation
49    use status
50 +  use interpolation
51   #ifdef IS_MPI
52    use mpiSimulation
53   #endif
# Line 64 | Line 65 | module electrostatic_module
65    !! these prefactors convert the multipole interactions into kcal / mol
66    !! all were computed assuming distances are measured in angstroms
67    !! Charge-Charge, assuming charges are measured in electrons
68 <  real(kind=dp), parameter :: pre11 = 332.0637778_dp
68 >  real(kind=dp), parameter :: pre11 = 332.0637778d0
69    !! Charge-Dipole, assuming charges are measured in electrons, and
70    !! dipoles are measured in debyes
71 <  real(kind=dp), parameter :: pre12 = 69.13373_dp
71 >  real(kind=dp), parameter :: pre12 = 69.13373d0
72    !! Dipole-Dipole, assuming dipoles are measured in debyes
73 <  real(kind=dp), parameter :: pre22 = 14.39325_dp
73 >  real(kind=dp), parameter :: pre22 = 14.39325d0
74    !! Charge-Quadrupole, assuming charges are measured in electrons, and
75    !! quadrupoles are measured in 10^-26 esu cm^2
76    !! This unit is also known affectionately as an esu centi-barn.
77 <  real(kind=dp), parameter :: pre14 = 69.13373_dp
77 >  real(kind=dp), parameter :: pre14 = 69.13373d0
78  
79    !! variables to handle different summation methods for long-range
80    !! electrostatics:
# Line 121 | Line 122 | module electrostatic_module
122    public :: setElectrostaticCutoffRadius
123    public :: setDampingAlpha
124    public :: setReactionFieldDielectric
125 +  public :: buildElectroSplines
126    public :: newElectrostaticType
127    public :: setCharge
128    public :: setDipoleMoment
# Line 133 | Line 135 | module electrostatic_module
135    public :: self_self
136    public :: rf_self_excludes
137  
138 +
139    type :: Electrostatic
140       integer :: c_ident
141       logical :: is_Charge = .false.
# Line 148 | Line 151 | contains
151  
152    type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap
153  
154 +  logical, save :: hasElectrostaticMap
155 +
156   contains
157  
158    subroutine setElectrostaticSummationMethod(the_ESM)
# Line 188 | Line 193 | contains
193      dielectric = thisDielectric
194      haveDielectric = .true.
195    end subroutine setReactionFieldDielectric
196 +
197 +  subroutine buildElectroSplines()
198 +  end subroutine buildElectroSplines
199  
200    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
201         is_SplitDipole, is_Quadrupole, is_Tap, status)
# Line 216 | Line 224 | contains
224            return
225         end if
226  
227 <       if (.not. allocated(ElectrostaticMap)) then
220 <          allocate(ElectrostaticMap(nAtypes))
221 <       endif
227 >       allocate(ElectrostaticMap(nAtypes))
228  
229      end if
230  
# Line 236 | Line 242 | contains
242      ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
243      ElectrostaticMap(myATID)%is_Tap = is_Tap
244  
245 +    hasElectrostaticMap = .true.
246 +
247    end subroutine newElectrostaticType
248  
249    subroutine setCharge(c_ident, charge, status)
# Line 247 | Line 255 | contains
255      status = 0
256      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
257  
258 <    if (.not.allocated(ElectrostaticMap)) then
258 >    if (.not.hasElectrostaticMap) then
259         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setCharge!")
260         status = -1
261         return
# Line 277 | Line 285 | contains
285      status = 0
286      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
287  
288 <    if (.not.allocated(ElectrostaticMap)) then
288 >    if (.not.hasElectrostaticMap) then
289         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setDipoleMoment!")
290         status = -1
291         return
# Line 307 | Line 315 | contains
315      status = 0
316      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
317  
318 <    if (.not.allocated(ElectrostaticMap)) then
318 >    if (.not.hasElectrostaticMap) then
319         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setSplitDipoleDistance!")
320         status = -1
321         return
# Line 337 | Line 345 | contains
345      status = 0
346      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
347  
348 <    if (.not.allocated(ElectrostaticMap)) then
348 >    if (.not.hasElectrostaticMap) then
349         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setQuadrupoleMoments!")
350         status = -1
351         return
# Line 368 | Line 376 | contains
376      integer :: localError
377      real(kind=dp) :: c
378  
379 <    if (.not.allocated(ElectrostaticMap)) then
379 >    if (.not.hasElectrostaticMap) then
380         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
381         return
382      end if
# Line 386 | Line 394 | contains
394      integer :: localError
395      real(kind=dp) :: dm
396  
397 <    if (.not.allocated(ElectrostaticMap)) then
397 >    if (.not.hasElectrostaticMap) then
398         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
399         return
400      end if
# Line 488 | Line 496 | contains
496      real (kind=dp) :: pot_term
497      real (kind=dp) :: preVal, rfVal
498      real (kind=dp) :: f13, f134
491 –
492 –    if (.not.allocated(ElectrostaticMap)) then
493 –       call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
494 –       return
495 –    end if
499  
500      if (.not.summationMethodChecked) then
501         call checkSummationMethod()
# Line 632 | Line 635 | contains
635         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
636      endif
637    
638 <    epot = 0.0_dp
639 <    dudx = 0.0_dp
640 <    dudy = 0.0_dp
641 <    dudz = 0.0_dp
638 >    epot = 0.0d0
639 >    dudx = 0.0d0
640 >    dudy = 0.0d0
641 >    dudz = 0.0d0
642  
643 <    dudux_i = 0.0_dp
644 <    duduy_i = 0.0_dp
645 <    duduz_i = 0.0_dp
643 >    dudux_i = 0.0d0
644 >    duduy_i = 0.0d0
645 >    duduz_i = 0.0d0
646  
647 <    dudux_j = 0.0_dp
648 <    duduy_j = 0.0_dp
649 <    duduz_j = 0.0_dp
647 >    dudux_j = 0.0d0
648 >    duduy_j = 0.0d0
649 >    duduz_j = 0.0d0
650  
651      if (i_is_Charge) then
652  
# Line 723 | Line 726 | contains
726  
727            else
728               if (j_is_SplitDipole) then
729 <                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
730 <                ri = 1.0_dp / BigR
729 >                BigR = sqrt(r2 + 0.25d0 * d_j * d_j)
730 >                ri = 1.0d0 / BigR
731                  scale = rij * ri
732               else
733                  ri = riji
734 <                scale = 1.0_dp
734 >                scale = 1.0d0
735               endif
736              
737               ri2 = ri * ri
# Line 775 | Line 778 | contains
778            cy2 = cy_j * cy_j
779            cz2 = cz_j * cz_j
780  
781 <          pref =  pre14 * q_i / 3.0_dp
782 <          pot_term = ri3*(qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
783 <               qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
784 <               qzz_j * (3.0_dp*cz2 - 1.0_dp))
781 >          pref =  pre14 * q_i / 3.0d0
782 >          pot_term = ri3*(qxx_j * (3.0d0*cx2 - 1.0d0) + &
783 >               qyy_j * (3.0d0*cy2 - 1.0d0) + &
784 >               qzz_j * (3.0d0*cz2 - 1.0d0))
785            vterm = pref * (pot_term*f1 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f2)
786            vpair = vpair + vterm
787            epot = epot + sw*vterm
788            
789            dudx = dudx - sw*pref*pot_term*riji*xhat*(5.0d0*f1 + f3) + &
790                 sw*pref*ri4 * ( &
791 <               qxx_j*(2.0_dp*cx_j*ux_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
792 <               qyy_j*(2.0_dp*cy_j*uy_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
793 <               qzz_j*(2.0_dp*cz_j*uz_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) ) &
791 >               qxx_j*(2.0d0*cx_j*ux_j(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) + &
792 >               qyy_j*(2.0d0*cy_j*uy_j(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) + &
793 >               qzz_j*(2.0d0*cz_j*uz_j(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) ) &
794                 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
795            dudy = dudy - sw*pref*pot_term*riji*yhat*(5.0d0*f1 + f3) + &
796                 sw*pref*ri4 * ( &
797 <               qxx_j*(2.0_dp*cx_j*ux_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
798 <               qyy_j*(2.0_dp*cy_j*uy_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
799 <               qzz_j*(2.0_dp*cz_j*uz_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) ) &
797 >               qxx_j*(2.0d0*cx_j*ux_j(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) + &
798 >               qyy_j*(2.0d0*cy_j*uy_j(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) + &
799 >               qzz_j*(2.0d0*cz_j*uz_j(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) ) &
800                 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
801            dudz = dudz - sw*pref*pot_term*riji*zhat*(5.0d0*f1 + f3) + &
802                 sw*pref*ri4 * ( &
803 <               qxx_j*(2.0_dp*cx_j*ux_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
804 <               qyy_j*(2.0_dp*cy_j*uy_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
805 <               qzz_j*(2.0_dp*cz_j*uz_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) ) &
803 >               qxx_j*(2.0d0*cx_j*ux_j(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) + &
804 >               qyy_j*(2.0d0*cy_j*uy_j(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) + &
805 >               qzz_j*(2.0d0*cz_j*uz_j(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) ) &
806                 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
807            
808 <          dudux_j(1) = dudux_j(1) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*xhat) &
808 >          dudux_j(1) = dudux_j(1) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*xhat) &
809                 * (3.0d0*f1 + f3) )
810 <          dudux_j(2) = dudux_j(2) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*yhat) &
810 >          dudux_j(2) = dudux_j(2) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*yhat) &
811                 * (3.0d0*f1 + f3) )
812 <          dudux_j(3) = dudux_j(3) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*zhat) &
812 >          dudux_j(3) = dudux_j(3) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*zhat) &
813                 * (3.0d0*f1 + f3) )
814            
815 <          duduy_j(1) = duduy_j(1) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*xhat) &
815 >          duduy_j(1) = duduy_j(1) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*xhat) &
816                 * (3.0d0*f1 + f3) )
817 <          duduy_j(2) = duduy_j(2) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*yhat) &
817 >          duduy_j(2) = duduy_j(2) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*yhat) &
818                 * (3.0d0*f1 + f3) )
819 <          duduy_j(3) = duduy_j(3) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*zhat) &
819 >          duduy_j(3) = duduy_j(3) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*zhat) &
820                 * (3.0d0*f1 + f3) )
821            
822 <          duduz_j(1) = duduz_j(1) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*xhat) &
822 >          duduz_j(1) = duduz_j(1) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*xhat) &
823                 * (3.0d0*f1 + f3) )
824 <          duduz_j(2) = duduz_j(2) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*yhat) &
824 >          duduz_j(2) = duduz_j(2) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*yhat) &
825                 * (3.0d0*f1 + f3) )
826 <          duduz_j(3) = duduz_j(3) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*zhat) &
826 >          duduz_j(3) = duduz_j(3) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*zhat) &
827                 * (3.0d0*f1 + f3) )
828            
829         endif
# Line 898 | Line 901 | contains
901  
902            else
903               if (i_is_SplitDipole) then
904 <                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
905 <                ri = 1.0_dp / BigR
904 >                BigR = sqrt(r2 + 0.25d0 * d_i * d_i)
905 >                ri = 1.0d0 / BigR
906                  scale = rij * ri
907               else
908                  ri = riji
909 <                scale = 1.0_dp
909 >                scale = 1.0d0
910               endif
911              
912               ri2 = ri * ri
# Line 977 | Line 980 | contains
980            else
981               if (i_is_SplitDipole) then
982                  if (j_is_SplitDipole) then
983 <                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
983 >                   BigR = sqrt(r2 + 0.25d0 * d_i * d_i + 0.25d0 * d_j * d_j)
984                  else
985 <                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
985 >                   BigR = sqrt(r2 + 0.25d0 * d_i * d_i)
986                  endif
987 <                ri = 1.0_dp / BigR
987 >                ri = 1.0d0 / BigR
988                  scale = rij * ri                
989               else
990                  if (j_is_SplitDipole) then
991 <                   BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
992 <                   ri = 1.0_dp / BigR
991 >                   BigR = sqrt(r2 + 0.25d0 * d_j * d_j)
992 >                   ri = 1.0d0 / BigR
993                     scale = rij * ri                            
994                  else                
995                     ri = riji
996 <                   scale = 1.0_dp
996 >                   scale = 1.0d0
997                  endif
998               endif
999              
# Line 1065 | Line 1068 | contains
1068            cy2 = cy_i * cy_i
1069            cz2 = cz_i * cz_i
1070  
1071 <          pref = pre14 * q_j / 3.0_dp
1072 <          pot_term = ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1073 <                            qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1074 <                            qzz_i * (3.0_dp*cz2 - 1.0_dp))
1071 >          pref = pre14 * q_j / 3.0d0
1072 >          pot_term = ri3 * (qxx_i * (3.0d0*cx2 - 1.0d0) + &
1073 >                            qyy_i * (3.0d0*cy2 - 1.0d0) + &
1074 >                            qzz_i * (3.0d0*cz2 - 1.0d0))
1075            vterm = pref * (pot_term*f1 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f2)
1076            vpair = vpair + vterm
1077            epot = epot + sw*vterm
1078            
1079            dudx = dudx - sw*pref*pot_term*riji*xhat*(5.0d0*f1 + f3) + &
1080                 sw*pref*ri4 * ( &
1081 <               qxx_i*(2.0_dp*cx_i*ux_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
1082 <               qyy_i*(2.0_dp*cy_i*uy_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
1083 <               qzz_i*(2.0_dp*cz_i*uz_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) ) &
1081 >               qxx_i*(2.0d0*cx_i*ux_i(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) + &
1082 >               qyy_i*(2.0d0*cy_i*uy_i(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) + &
1083 >               qzz_i*(2.0d0*cz_i*uz_i(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) ) &
1084                 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1085            dudy = dudy - sw*pref*pot_term*riji*yhat*(5.0d0*f1 + f3) + &
1086                 sw*pref*ri4 * ( &
1087 <               qxx_i*(2.0_dp*cx_i*ux_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
1088 <               qyy_i*(2.0_dp*cy_i*uy_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
1089 <               qzz_i*(2.0_dp*cz_i*uz_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) ) &
1087 >               qxx_i*(2.0d0*cx_i*ux_i(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) + &
1088 >               qyy_i*(2.0d0*cy_i*uy_i(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) + &
1089 >               qzz_i*(2.0d0*cz_i*uz_i(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) ) &
1090                 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1091            dudz = dudz - sw*pref*pot_term*riji*zhat*(5.0d0*f1 + f3) + &
1092                 sw*pref*ri4 * ( &
1093 <               qxx_i*(2.0_dp*cx_i*ux_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
1094 <               qyy_i*(2.0_dp*cy_i*uy_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
1095 <               qzz_i*(2.0_dp*cz_i*uz_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) ) &
1093 >               qxx_i*(2.0d0*cx_i*ux_i(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) + &
1094 >               qyy_i*(2.0d0*cy_i*uy_i(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) + &
1095 >               qzz_i*(2.0d0*cz_i*uz_i(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) ) &
1096                 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1097            
1098 <          dudux_i(1) = dudux_i(1) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*xhat) &
1098 >          dudux_i(1) = dudux_i(1) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*xhat) &
1099                 * (3.0d0*f1 + f3) )
1100 <          dudux_i(2) = dudux_i(2) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*yhat) &
1100 >          dudux_i(2) = dudux_i(2) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*yhat) &
1101                 * (3.0d0*f1 + f3) )
1102 <          dudux_i(3) = dudux_i(3) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*zhat) &
1102 >          dudux_i(3) = dudux_i(3) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*zhat) &
1103                 * (3.0d0*f1 + f3) )
1104            
1105 <          duduy_i(1) = duduy_i(1) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*xhat) &
1105 >          duduy_i(1) = duduy_i(1) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*xhat) &
1106                 * (3.0d0*f1 + f3) )
1107 <          duduy_i(2) = duduy_i(2) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*yhat) &
1107 >          duduy_i(2) = duduy_i(2) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*yhat) &
1108                 * (3.0d0*f1 + f3) )
1109 <          duduy_i(3) = duduy_i(3) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*zhat) &
1109 >          duduy_i(3) = duduy_i(3) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*zhat) &
1110                 * (3.0d0*f1 + f3) )
1111            
1112 <          duduz_i(1) = duduz_i(1) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*xhat) &
1112 >          duduz_i(1) = duduz_i(1) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*xhat) &
1113                 * (3.0d0*f1 + f3) )
1114 <          duduz_i(2) = duduz_i(2) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*yhat) &
1114 >          duduz_i(2) = duduz_i(2) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*yhat) &
1115                 * (3.0d0*f1 + f3) )
1116 <          duduz_i(3) = duduz_i(3) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*zhat) &
1116 >          duduz_i(3) = duduz_i(3) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*zhat) &
1117                 * (3.0d0*f1 + f3) )
1118  
1119         endif
# Line 1279 | Line 1282 | contains
1282            c1 = getCharge(atid1)
1283            
1284            if (screeningMethod .eq. DAMPED) then
1285 <             mypot = mypot - (f0c * rcuti * 0.5_dp + &
1285 >             mypot = mypot - (f0c * rcuti * 0.5d0 + &
1286                    dampingAlpha*invRootPi) * c1 * c1    
1287              
1288            else            
1289 <             mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1289 >             mypot = mypot - (rcuti * 0.5d0 * c1 * c1)
1290              
1291            endif
1292         endif

Diff Legend

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