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 2717 by gezelter, Mon Apr 17 21:49:12 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    !! variables to handle different summation methods for long-range
80    !! electrostatics:
# Line 151 | 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 222 | Line 224 | contains
224            return
225         end if
226  
227 <       if (.not. allocated(ElectrostaticMap)) then
226 <          allocate(ElectrostaticMap(nAtypes))
227 <       endif
227 >       allocate(ElectrostaticMap(nAtypes))
228  
229      end if
230  
# Line 242 | 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 253 | 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 283 | 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 313 | 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 343 | 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 374 | 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 392 | 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 494 | Line 496 | contains
496      real (kind=dp) :: pot_term
497      real (kind=dp) :: preVal, rfVal
498      real (kind=dp) :: f13, f134
497 –
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
499  
500      if (.not.summationMethodChecked) then
501         call checkSummationMethod()
# Line 638 | 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 729 | 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 781 | 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 904 | 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 983 | 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 1071 | 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 1285 | 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