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 2715 by chrisfen, Sun Apr 16 02:51:16 2006 UTC vs.
Revision 2722 by gezelter, Thu Apr 20 18:24:24 2006 UTC

# Line 65 | 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 +  real(kind=dp), parameter :: zero = 0.0d0
80 +
81    !! variables to handle different summation methods for long-range
82    !! electrostatics:
83    integer, save :: summationMethod = NONE
# Line 151 | Line 153 | contains
153  
154    type(Electrostatic), dimension(:), allocatable :: ElectrostaticMap
155  
156 +  logical, save :: hasElectrostaticMap
157 +
158   contains
159  
160    subroutine setElectrostaticSummationMethod(the_ESM)
# Line 222 | Line 226 | contains
226            return
227         end if
228  
229 <       if (.not. allocated(ElectrostaticMap)) then
226 <          allocate(ElectrostaticMap(nAtypes))
227 <       endif
229 >       allocate(ElectrostaticMap(nAtypes))
230  
231      end if
232  
# Line 242 | Line 244 | contains
244      ElectrostaticMap(myATID)%is_Quadrupole = is_Quadrupole
245      ElectrostaticMap(myATID)%is_Tap = is_Tap
246  
247 +    hasElectrostaticMap = .true.
248 +
249    end subroutine newElectrostaticType
250  
251    subroutine setCharge(c_ident, charge, status)
# Line 253 | Line 257 | contains
257      status = 0
258      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
259  
260 <    if (.not.allocated(ElectrostaticMap)) then
260 >    if (.not.hasElectrostaticMap) then
261         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setCharge!")
262         status = -1
263         return
# Line 283 | Line 287 | contains
287      status = 0
288      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
289  
290 <    if (.not.allocated(ElectrostaticMap)) then
290 >    if (.not.hasElectrostaticMap) then
291         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setDipoleMoment!")
292         status = -1
293         return
# Line 313 | Line 317 | contains
317      status = 0
318      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
319  
320 <    if (.not.allocated(ElectrostaticMap)) then
320 >    if (.not.hasElectrostaticMap) then
321         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setSplitDipoleDistance!")
322         status = -1
323         return
# Line 343 | Line 347 | contains
347      status = 0
348      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
349  
350 <    if (.not.allocated(ElectrostaticMap)) then
350 >    if (.not.hasElectrostaticMap) then
351         call handleError("electrostatic", "no ElectrostaticMap was present before first call of setQuadrupoleMoments!")
352         status = -1
353         return
# Line 374 | Line 378 | contains
378      integer :: localError
379      real(kind=dp) :: c
380  
381 <    if (.not.allocated(ElectrostaticMap)) then
381 >    if (.not.hasElectrostaticMap) then
382         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getCharge!")
383         return
384      end if
# Line 392 | Line 396 | contains
396      integer :: localError
397      real(kind=dp) :: dm
398  
399 <    if (.not.allocated(ElectrostaticMap)) then
399 >    if (.not.hasElectrostaticMap) then
400         call handleError("electrostatic", "no ElectrostaticMap was present before first call of getDipoleMoment!")
401         return
402      end if
# Line 495 | Line 499 | contains
499      real (kind=dp) :: preVal, rfVal
500      real (kind=dp) :: f13, f134
501  
498 –    if (.not.allocated(ElectrostaticMap)) then
499 –       call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
500 –       return
501 –    end if
502 –
502      if (.not.summationMethodChecked) then
503         call checkSummationMethod()
504      endif
# Line 553 | Line 552 | contains
552         if (i_is_SplitDipole) then
553            d_i = ElectrostaticMap(me1)%split_dipole_distance
554         endif
555 <
555 >       duduz_i = zero
556      endif
557  
558      if (i_is_Quadrupole) then
# Line 584 | Line 583 | contains
583         cx_i = ux_i(1)*xhat + ux_i(2)*yhat + ux_i(3)*zhat
584         cy_i = uy_i(1)*xhat + uy_i(2)*yhat + uy_i(3)*zhat
585         cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
586 +       dudux_i = zero
587 +       duduy_i = zero
588 +       duduz_i = zero
589      endif
590  
591      if (j_is_Charge) then
# Line 606 | Line 608 | contains
608         if (j_is_SplitDipole) then
609            d_j = ElectrostaticMap(me2)%split_dipole_distance
610         endif
611 +       duduz_j = zero
612      endif
613  
614      if (j_is_Quadrupole) then
# Line 636 | Line 639 | contains
639         cx_j = ux_j(1)*xhat + ux_j(2)*yhat + ux_j(3)*zhat
640         cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
641         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
642 +       dudux_j = zero
643 +       duduy_j = zero
644 +       duduz_j = zero
645      endif
646    
647 <    epot = 0.0_dp
648 <    dudx = 0.0_dp
649 <    dudy = 0.0_dp
650 <    dudz = 0.0_dp
645 <
646 <    dudux_i = 0.0_dp
647 <    duduy_i = 0.0_dp
648 <    duduz_i = 0.0_dp
649 <
650 <    dudux_j = 0.0_dp
651 <    duduy_j = 0.0_dp
652 <    duduz_j = 0.0_dp
647 >    epot = zero
648 >    dudx = zero
649 >    dudy = zero
650 >    dudz = zero  
651  
652      if (i_is_Charge) then
653  
# Line 729 | Line 727 | contains
727  
728            else
729               if (j_is_SplitDipole) then
730 <                BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
731 <                ri = 1.0_dp / BigR
730 >                BigR = sqrt(r2 + 0.25d0 * d_j * d_j)
731 >                ri = 1.0d0 / BigR
732                  scale = rij * ri
733               else
734                  ri = riji
735 <                scale = 1.0_dp
735 >                scale = 1.0d0
736               endif
737              
738               ri2 = ri * ri
# Line 781 | Line 779 | contains
779            cy2 = cy_j * cy_j
780            cz2 = cz_j * cz_j
781  
782 <          pref =  pre14 * q_i / 3.0_dp
783 <          pot_term = ri3*(qxx_j * (3.0_dp*cx2 - 1.0_dp) + &
784 <               qyy_j * (3.0_dp*cy2 - 1.0_dp) + &
785 <               qzz_j * (3.0_dp*cz2 - 1.0_dp))
782 >          pref =  pre14 * q_i / 3.0d0
783 >          pot_term = ri3*(qxx_j * (3.0d0*cx2 - 1.0d0) + &
784 >               qyy_j * (3.0d0*cy2 - 1.0d0) + &
785 >               qzz_j * (3.0d0*cz2 - 1.0d0))
786            vterm = pref * (pot_term*f1 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f2)
787            vpair = vpair + vterm
788            epot = epot + sw*vterm
789            
790            dudx = dudx - sw*pref*pot_term*riji*xhat*(5.0d0*f1 + f3) + &
791                 sw*pref*ri4 * ( &
792 <               qxx_j*(2.0_dp*cx_j*ux_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
793 <               qyy_j*(2.0_dp*cy_j*uy_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
794 <               qzz_j*(2.0_dp*cz_j*uz_j(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) ) &
792 >               qxx_j*(2.0d0*cx_j*ux_j(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) + &
793 >               qyy_j*(2.0d0*cy_j*uy_j(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) + &
794 >               qzz_j*(2.0d0*cz_j*uz_j(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) ) &
795                 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
796            dudy = dudy - sw*pref*pot_term*riji*yhat*(5.0d0*f1 + f3) + &
797                 sw*pref*ri4 * ( &
798 <               qxx_j*(2.0_dp*cx_j*ux_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
799 <               qyy_j*(2.0_dp*cy_j*uy_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
800 <               qzz_j*(2.0_dp*cz_j*uz_j(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) ) &
798 >               qxx_j*(2.0d0*cx_j*ux_j(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) + &
799 >               qyy_j*(2.0d0*cy_j*uy_j(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) + &
800 >               qzz_j*(2.0d0*cz_j*uz_j(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) ) &
801                 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
802            dudz = dudz - sw*pref*pot_term*riji*zhat*(5.0d0*f1 + f3) + &
803                 sw*pref*ri4 * ( &
804 <               qxx_j*(2.0_dp*cx_j*ux_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
805 <               qyy_j*(2.0_dp*cy_j*uy_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
806 <               qzz_j*(2.0_dp*cz_j*uz_j(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) ) &
804 >               qxx_j*(2.0d0*cx_j*ux_j(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) + &
805 >               qyy_j*(2.0d0*cy_j*uy_j(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) + &
806 >               qzz_j*(2.0d0*cz_j*uz_j(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) ) &
807                 + (qxx_j*cx2 + qyy_j*cy2 + qzz_j*cz2)*f4
808            
809 <          dudux_j(1) = dudux_j(1) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*xhat) &
809 >          dudux_j(1) = dudux_j(1) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*xhat) &
810                 * (3.0d0*f1 + f3) )
811 <          dudux_j(2) = dudux_j(2) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*yhat) &
811 >          dudux_j(2) = dudux_j(2) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*yhat) &
812                 * (3.0d0*f1 + f3) )
813 <          dudux_j(3) = dudux_j(3) + sw*pref*ri3*( (qxx_j*2.0_dp*cx_j*zhat) &
813 >          dudux_j(3) = dudux_j(3) + sw*pref*ri3*( (qxx_j*2.0d0*cx_j*zhat) &
814                 * (3.0d0*f1 + f3) )
815            
816 <          duduy_j(1) = duduy_j(1) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*xhat) &
816 >          duduy_j(1) = duduy_j(1) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*xhat) &
817                 * (3.0d0*f1 + f3) )
818 <          duduy_j(2) = duduy_j(2) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*yhat) &
818 >          duduy_j(2) = duduy_j(2) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*yhat) &
819                 * (3.0d0*f1 + f3) )
820 <          duduy_j(3) = duduy_j(3) + sw*pref*ri3*( (qyy_j*2.0_dp*cy_j*zhat) &
820 >          duduy_j(3) = duduy_j(3) + sw*pref*ri3*( (qyy_j*2.0d0*cy_j*zhat) &
821                 * (3.0d0*f1 + f3) )
822            
823 <          duduz_j(1) = duduz_j(1) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*xhat) &
823 >          duduz_j(1) = duduz_j(1) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*xhat) &
824                 * (3.0d0*f1 + f3) )
825 <          duduz_j(2) = duduz_j(2) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*yhat) &
825 >          duduz_j(2) = duduz_j(2) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*yhat) &
826                 * (3.0d0*f1 + f3) )
827 <          duduz_j(3) = duduz_j(3) + sw*pref*ri3*( (qzz_j*2.0_dp*cz_j*zhat) &
827 >          duduz_j(3) = duduz_j(3) + sw*pref*ri3*( (qzz_j*2.0d0*cz_j*zhat) &
828                 * (3.0d0*f1 + f3) )
829            
830         endif
# Line 904 | Line 902 | contains
902  
903            else
904               if (i_is_SplitDipole) then
905 <                BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
906 <                ri = 1.0_dp / BigR
905 >                BigR = sqrt(r2 + 0.25d0 * d_i * d_i)
906 >                ri = 1.0d0 / BigR
907                  scale = rij * ri
908               else
909                  ri = riji
910 <                scale = 1.0_dp
910 >                scale = 1.0d0
911               endif
912              
913               ri2 = ri * ri
# Line 983 | Line 981 | contains
981            else
982               if (i_is_SplitDipole) then
983                  if (j_is_SplitDipole) then
984 <                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i + 0.25_dp * d_j * d_j)
984 >                   BigR = sqrt(r2 + 0.25d0 * d_i * d_i + 0.25d0 * d_j * d_j)
985                  else
986 <                   BigR = sqrt(r2 + 0.25_dp * d_i * d_i)
986 >                   BigR = sqrt(r2 + 0.25d0 * d_i * d_i)
987                  endif
988 <                ri = 1.0_dp / BigR
988 >                ri = 1.0d0 / BigR
989                  scale = rij * ri                
990               else
991                  if (j_is_SplitDipole) then
992 <                   BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
993 <                   ri = 1.0_dp / BigR
992 >                   BigR = sqrt(r2 + 0.25d0 * d_j * d_j)
993 >                   ri = 1.0d0 / BigR
994                     scale = rij * ri                            
995                  else                
996                     ri = riji
997 <                   scale = 1.0_dp
997 >                   scale = 1.0d0
998                  endif
999               endif
1000              
# Line 1071 | Line 1069 | contains
1069            cy2 = cy_i * cy_i
1070            cz2 = cz_i * cz_i
1071  
1072 <          pref = pre14 * q_j / 3.0_dp
1073 <          pot_term = ri3 * (qxx_i * (3.0_dp*cx2 - 1.0_dp) + &
1074 <                            qyy_i * (3.0_dp*cy2 - 1.0_dp) + &
1075 <                            qzz_i * (3.0_dp*cz2 - 1.0_dp))
1072 >          pref = pre14 * q_j / 3.0d0
1073 >          pot_term = ri3 * (qxx_i * (3.0d0*cx2 - 1.0d0) + &
1074 >                            qyy_i * (3.0d0*cy2 - 1.0d0) + &
1075 >                            qzz_i * (3.0d0*cz2 - 1.0d0))
1076            vterm = pref * (pot_term*f1 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f2)
1077            vpair = vpair + vterm
1078            epot = epot + sw*vterm
1079            
1080            dudx = dudx - sw*pref*pot_term*riji*xhat*(5.0d0*f1 + f3) + &
1081                 sw*pref*ri4 * ( &
1082 <               qxx_i*(2.0_dp*cx_i*ux_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
1083 <               qyy_i*(2.0_dp*cy_i*uy_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) + &
1084 <               qzz_i*(2.0_dp*cz_i*uz_i(1)*(3.0d0*f1 + f3) - 2.0_dp*xhat*f1) ) &
1082 >               qxx_i*(2.0d0*cx_i*ux_i(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) + &
1083 >               qyy_i*(2.0d0*cy_i*uy_i(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) + &
1084 >               qzz_i*(2.0d0*cz_i*uz_i(1)*(3.0d0*f1 + f3) - 2.0d0*xhat*f1) ) &
1085                 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1086            dudy = dudy - sw*pref*pot_term*riji*yhat*(5.0d0*f1 + f3) + &
1087                 sw*pref*ri4 * ( &
1088 <               qxx_i*(2.0_dp*cx_i*ux_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
1089 <               qyy_i*(2.0_dp*cy_i*uy_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) + &
1090 <               qzz_i*(2.0_dp*cz_i*uz_i(2)*(3.0d0*f1 + f3) - 2.0_dp*yhat*f1) ) &
1088 >               qxx_i*(2.0d0*cx_i*ux_i(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) + &
1089 >               qyy_i*(2.0d0*cy_i*uy_i(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) + &
1090 >               qzz_i*(2.0d0*cz_i*uz_i(2)*(3.0d0*f1 + f3) - 2.0d0*yhat*f1) ) &
1091                 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1092            dudz = dudz - sw*pref*pot_term*riji*zhat*(5.0d0*f1 + f3) + &
1093                 sw*pref*ri4 * ( &
1094 <               qxx_i*(2.0_dp*cx_i*ux_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
1095 <               qyy_i*(2.0_dp*cy_i*uy_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) + &
1096 <               qzz_i*(2.0_dp*cz_i*uz_i(3)*(3.0d0*f1 + f3) - 2.0_dp*zhat*f1) ) &
1094 >               qxx_i*(2.0d0*cx_i*ux_i(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) + &
1095 >               qyy_i*(2.0d0*cy_i*uy_i(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) + &
1096 >               qzz_i*(2.0d0*cz_i*uz_i(3)*(3.0d0*f1 + f3) - 2.0d0*zhat*f1) ) &
1097                 + (qxx_i*cx2 + qyy_i*cy2 + qzz_i*cz2)*f4
1098            
1099 <          dudux_i(1) = dudux_i(1) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*xhat) &
1099 >          dudux_i(1) = dudux_i(1) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*xhat) &
1100                 * (3.0d0*f1 + f3) )
1101 <          dudux_i(2) = dudux_i(2) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*yhat) &
1101 >          dudux_i(2) = dudux_i(2) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*yhat) &
1102                 * (3.0d0*f1 + f3) )
1103 <          dudux_i(3) = dudux_i(3) + sw*pref*( ri3*(qxx_i*2.0_dp*cx_i*zhat) &
1103 >          dudux_i(3) = dudux_i(3) + sw*pref*( ri3*(qxx_i*2.0d0*cx_i*zhat) &
1104                 * (3.0d0*f1 + f3) )
1105            
1106 <          duduy_i(1) = duduy_i(1) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*xhat) &
1106 >          duduy_i(1) = duduy_i(1) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*xhat) &
1107                 * (3.0d0*f1 + f3) )
1108 <          duduy_i(2) = duduy_i(2) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*yhat) &
1108 >          duduy_i(2) = duduy_i(2) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*yhat) &
1109                 * (3.0d0*f1 + f3) )
1110 <          duduy_i(3) = duduy_i(3) + sw*pref*( ri3*(qyy_i*2.0_dp*cy_i*zhat) &
1110 >          duduy_i(3) = duduy_i(3) + sw*pref*( ri3*(qyy_i*2.0d0*cy_i*zhat) &
1111                 * (3.0d0*f1 + f3) )
1112            
1113 <          duduz_i(1) = duduz_i(1) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*xhat) &
1116 <               * (3.0d0*f1 + f3) )
1117 <          duduz_i(2) = duduz_i(2) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*yhat) &
1113 >          duduz_i(1) = duduz_i(1) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*xhat) &
1114                 * (3.0d0*f1 + f3) )
1115 <          duduz_i(3) = duduz_i(3) + sw*pref*( ri3*(qzz_i*2.0_dp*cz_i*zhat) &
1115 >          duduz_i(2) = duduz_i(2) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*yhat) &
1116                 * (3.0d0*f1 + f3) )
1117 +          duduz_i(3) = duduz_i(3) + sw*pref*( ri3*(qzz_i*2.0d0*cz_i*zhat) &
1118 +               * (3.0d0*f1 + f3) )
1119  
1120         endif
1121      endif
# Line 1285 | Line 1283 | contains
1283            c1 = getCharge(atid1)
1284            
1285            if (screeningMethod .eq. DAMPED) then
1286 <             mypot = mypot - (f0c * rcuti * 0.5_dp + &
1286 >             mypot = mypot - (f0c * rcuti * 0.5d0 + &
1287                    dampingAlpha*invRootPi) * c1 * c1    
1288              
1289            else            
1290 <             mypot = mypot - (rcuti * 0.5_dp * c1 * c1)
1290 >             mypot = mypot - (rcuti * 0.5d0 * c1 * c1)
1291              
1292            endif
1293         endif
# Line 1329 | Line 1327 | contains
1327         call checkSummationMethod()
1328      endif
1329  
1330 <    dudx = 0.0d0
1331 <    dudy = 0.0d0
1332 <    dudz = 0.0d0
1330 >    dudx = zero
1331 >    dudy = zero
1332 >    dudz = zero
1333  
1334      riji = 1.0d0/rij
1335  

Diff Legend

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