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 2339 by chrisfen, Wed Sep 28 18:47:17 2005 UTC vs.
Revision 2395 by chrisfen, Mon Oct 24 14:06:36 2005 UTC

# Line 54 | Line 54 | module electrostatic_module
54  
55    PRIVATE
56  
57 +
58   #define __FORTRAN90
59 + #include "UseTheForce/DarkSide/fInteractionMap.h"
60   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
61  
62 +
63    !! these prefactors convert the multipole interactions into kcal / mol
64    !! all were computed assuming distances are measured in angstroms
65    !! Charge-Charge, assuming charges are measured in electrons
# Line 75 | Line 78 | module electrostatic_module
78    integer, save :: summationMethod = NONE
79    logical, save :: summationMethodChecked = .false.
80    real(kind=DP), save :: defaultCutoff = 0.0_DP
81 +  real(kind=DP), save :: defaultCutoff2 = 0.0_DP
82    logical, save :: haveDefaultCutoff = .false.
83    real(kind=DP), save :: dampingAlpha = 0.0_DP
84    logical, save :: haveDampingAlpha = .false.
85 <  real(kind=DP), save :: dielectric = 0.0_DP
85 >  real(kind=DP), save :: dielectric = 1.0_DP
86    logical, save :: haveDielectric = .false.
87    real(kind=DP), save :: constERFC = 0.0_DP
88    real(kind=DP), save :: constEXP = 0.0_DP
89    logical, save :: haveDWAconstants = .false.
90 <  real(kind=dp), save :: rcuti = 0.0_dp
91 <  real(kind=dp), save :: rcuti2 = 0.0_dp
92 <  real(kind=dp), save :: rcuti3 = 0.0_dp
93 <  real(kind=dp), save :: rcuti4 = 0.0_dp
94 <  real(kind=dp), save :: alphaPi = 0.0_dp
95 <  real(kind=dp), save :: invRootPi = 0.0_dp
96 <  
90 >  real(kind=dp), save :: rcuti = 0.0_DP
91 >  real(kind=dp), save :: rcuti2 = 0.0_DP
92 >  real(kind=dp), save :: rcuti3 = 0.0_DP
93 >  real(kind=dp), save :: rcuti4 = 0.0_DP
94 >  real(kind=dp), save :: alphaPi = 0.0_DP
95 >  real(kind=dp), save :: invRootPi = 0.0_DP
96 >  real(kind=dp), save :: rrf = 1.0_DP
97 >  real(kind=dp), save :: rt = 1.0_DP
98 >  real(kind=dp), save :: rrfsq = 1.0_DP
99 >  real(kind=dp), save :: preRF = 0.0_DP
100 >  real(kind=dp), save :: preRF2 = 0.0_DP
101 >  logical, save :: preRFCalculated = .false.
102 >
103   #ifdef __IFC
104   ! error function for ifc version > 7.
105    double precision, external :: derfc
# Line 99 | Line 109 | module electrostatic_module
109    public :: setElectrostaticCutoffRadius
110    public :: setDampedWolfAlpha
111    public :: setReactionFieldDielectric
112 +  public :: setReactionFieldPrefactor
113    public :: newElectrostaticType
114    public :: setCharge
115    public :: setDipoleMoment
# Line 107 | Line 118 | module electrostatic_module
118    public :: doElectrostaticPair
119    public :: getCharge
120    public :: getDipoleMoment
110 –  public :: pre22
121    public :: destroyElectrostaticTypes
122 +  public :: rf_self_self
123  
124    type :: Electrostatic
125       integer :: c_ident
# Line 138 | Line 149 | contains
149  
150    end subroutine setElectrostaticSummationMethod
151  
152 <  subroutine setElectrostaticCutoffRadius(thisRcut)
152 >  subroutine setElectrostaticCutoffRadius(thisRcut, thisRsw)
153      real(kind=dp), intent(in) :: thisRcut
154 +    real(kind=dp), intent(in) :: thisRsw
155      defaultCutoff = thisRcut
156 +    rrf = defaultCutoff
157 +    rt = thisRsw
158      haveDefaultCutoff = .true.
159    end subroutine setElectrostaticCutoffRadius
160  
# Line 155 | Line 169 | contains
169      dielectric = thisDielectric
170      haveDielectric = .true.
171    end subroutine setReactionFieldDielectric
172 +
173 +  subroutine setReactionFieldPrefactor
174 +    if (haveDefaultCutoff .and. haveDielectric) then
175 +       defaultCutoff2 = defaultCutoff*defaultCutoff
176 +       preRF = (dielectric-1.0d0) / &
177 +            ((2.0d0*dielectric+1.0d0)*defaultCutoff2*defaultCutoff)
178 +       preRF2 = 2.0d0*preRF
179 +       preRFCalculated = .true.
180 +    else if (.not.haveDefaultCutoff) then
181 +       call handleError("setReactionFieldPrefactor", "Default cutoff not set")
182 +    else
183 +       call handleError("setReactionFieldPrefactor", "Dielectric not set")
184 +    endif
185 +  end subroutine setReactionFieldPrefactor
186  
187    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
188         is_SplitDipole, is_Quadrupole, is_Tap, status)
# Line 392 | Line 420 | contains
420            constERFC = derfc(dampingAlpha*defaultCutoff)
421            invRootPi = 0.56418958354775628695d0
422            alphaPi = 2*dampingAlpha*invRootPi
423 <          
423 >  
424            haveDWAconstants = .true.
425         endif
426      endif
# Line 448 | Line 476 | contains
476      real (kind=dp) :: dudx, dudy, dudz
477      real (kind=dp) :: scale, sc2, bigR
478      real (kind=dp) :: varERFC, varEXP
479 +    real (kind=dp) :: limScale
480 +    real (kind=dp) :: preVal, rfVal
481  
482      if (.not.allocated(ElectrostaticMap)) then
483         call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!")
# Line 456 | Line 486 | contains
486  
487      if (.not.summationMethodChecked) then
488         call checkSummationMethod()
459 –      
489      endif
490  
491 +    if (.not.preRFCalculated) then
492 +       call setReactionFieldPrefactor()
493 +    endif
494  
495   #ifdef IS_MPI
496      me1 = atid_Row(atom1)
# Line 471 | Line 503 | contains
503      !! some variables we'll need independent of electrostatic type:
504  
505      riji = 1.0d0 / rij
506 <
506 >  
507      xhat = d(1) * riji
508      yhat = d(2) * riji
509      zhat = d(3) * riji
# Line 612 | Line 644 | contains
644         if (j_is_Charge) then
645  
646            if (summationMethod .eq. UNDAMPED_WOLF) then
615 –
647               vterm = pre11 * q_i * q_j * (riji - rcuti)
648               vpair = vpair + vterm
649               epot = epot + sw*vterm
650              
651 <             dudr  = -sw*pre11*q_i*q_j * (riji*riji*riji - rcuti2*rcuti)
651 >             dudr  = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji
652              
653               dudx = dudx + dudr * d(1)
654               dudy = dudy + dudr * d(2)
655               dudz = dudz + dudr * d(3)
656  
657            elseif (summationMethod .eq. DAMPED_WOLF) then
627 –
658               varERFC = derfc(dampingAlpha*rij)
659               varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij)
660               vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti)
661               vpair = vpair + vterm
662               epot = epot + sw*vterm
663              
664 <             dudr  = -sw*pre11*q_i*q_j * ( riji*(varERFC*riji*riji &
665 <                                                 + alphaPi*varEXP) &
666 <                                         - rcuti*(constERFC*rcuti2 &
667 <                                                 + alphaPi*constEXP) )
664 >             dudr  = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji &
665 >                                                  + alphaPi*varEXP) &
666 >                                                 - (constERFC*rcuti2 &
667 >                                                    + alphaPi*constEXP)) )
668              
669               dudx = dudx + dudr * d(1)
670               dudy = dudy + dudr * d(2)
671               dudz = dudz + dudr * d(3)
672  
673 <          else
673 >          elseif (summationMethod .eq. REACTION_FIELD) then
674 >             preVal = pre11 * q_i * q_j
675 >             rfVal = preRF*rij*rij
676 >             vterm = preVal * ( riji + rfVal )
677 >             vpair = vpair + vterm
678 >             epot = epot + sw*vterm
679 >            
680 >             dudr  = sw * preVal * ( 2.0d0*rfVal - riji )*riji
681 >            
682 >             dudx = dudx + dudr * xhat
683 >             dudy = dudy + dudr * yhat
684 >             dudz = dudz + dudr * zhat
685  
686 +          else
687               vterm = pre11 * q_i * q_j * riji
688               vpair = vpair + vterm
689               epot = epot + sw*vterm
# Line 684 | Line 726 | contains
726               duduz_j(2) = duduz_j(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
727               duduz_j(3) = duduz_j(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
728  
729 +          elseif (summationMethod .eq. REACTION_FIELD) then
730 +             ri2 = ri * ri
731 +             ri3 = ri2 * ri
732 +    
733 +             pref = pre12 * q_i * mu_j
734 +             vterm = - pref * ct_j * ( ri2 - preRF2*rij )
735 +             vpair = vpair + vterm
736 +             epot = epot + sw*vterm
737 +            
738 +             !! this has a + sign in the () because the rij vector is
739 +             !! r_j - r_i and the charge-dipole potential takes the origin
740 +             !! as the point dipole, which is atom j in this case.
741 +            
742 +             dudx = dudx - sw*pref*( ri3*(uz_j(1) - 3.0d0*ct_j*xhat) - &
743 +                                     preRF2*uz_j(1) )
744 +             dudy = dudy - sw*pref*( ri3*(uz_j(2) - 3.0d0*ct_j*yhat) - &
745 +                                     preRF2*uz_j(2) )
746 +             dudz = dudz - sw*pref*( ri3*(uz_j(3) - 3.0d0*ct_j*zhat) - &
747 +                                     preRF2*uz_j(3) )        
748 +             duduz_j(1) = duduz_j(1) - sw*pref * xhat * ( ri2 - preRF2*rij )
749 +             duduz_j(2) = duduz_j(2) - sw*pref * yhat * ( ri2 - preRF2*rij )
750 +             duduz_j(3) = duduz_j(3) - sw*pref * zhat * ( ri2 - preRF2*rij )
751 +
752            else
753               if (j_is_SplitDipole) then
754                  BigR = sqrt(r2 + 0.25_dp * d_j * d_j)
# Line 849 | Line 914 | contains
914               duduz_i(1) = duduz_i(1) - sw*pref*( ri2*xhat - d(1)*rcuti3 )
915               duduz_i(2) = duduz_i(2) - sw*pref*( ri2*yhat - d(2)*rcuti3 )
916               duduz_i(3) = duduz_i(3) - sw*pref*( ri2*zhat - d(3)*rcuti3 )
917 +
918 +          elseif (summationMethod .eq. REACTION_FIELD) then
919 +             ri2 = ri * ri
920 +             ri3 = ri2 * ri
921 +
922 +             pref = pre12 * q_j * mu_i
923 +             vterm = pref * ct_i * ( ri2 - preRF*rij )
924 +             vpair = vpair + vterm
925 +             epot = epot + sw*vterm
926 +            
927 +             dudx = dudx + sw*pref * ri3 * ( uz_i(1) - 3.0d0*ct_i*xhat - &
928 +                                             preRF*uz_i(1) )
929 +             dudy = dudy + sw*pref * ri3 * ( uz_i(2) - 3.0d0*ct_i*yhat - &
930 +                                             preRF*uz_i(2) )
931 +             dudz = dudz + sw*pref * ri3 * ( uz_i(3) - 3.0d0*ct_i*zhat - &
932 +                                             preRF*uz_i(3) )
933 +            
934 +             duduz_i(1) = duduz_i(1) + sw*pref * xhat * ( ri2 - preRF*rij )
935 +             duduz_i(2) = duduz_i(2) + sw*pref * yhat * ( ri2 - preRF*rij )
936 +             duduz_i(3) = duduz_i(3) + sw*pref * zhat * ( ri2 - preRF*rij )
937  
938            else
939               if (i_is_SplitDipole) then
# Line 918 | Line 1003 | contains
1003                    - rcuti3*(uz_i(2) - 3.0d0*ct_i*d(2)*rcuti))
1004               duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1005                    - rcuti3*(uz_i(3) - 3.0d0*ct_i*d(3)*rcuti))
1006 +
1007 +         elseif (summationMethod .eq. REACTION_FIELD) then
1008 +             ct_ij = uz_i(1)*uz_j(1) + uz_i(2)*uz_j(2) + uz_i(3)*uz_j(3)
1009 +
1010 +             ri2 = riji * riji
1011 +             ri3 = ri2 * riji
1012 +             ri4 = ri2 * ri2
1013 +
1014 +             pref = pre22 * mu_i * mu_j
1015 +              
1016 +             vterm = pref*( ri3*(ct_ij - 3.0d0 * ct_i * ct_j) - &
1017 +                  preRF2*ct_ij )
1018 +             vpair = vpair + vterm
1019 +             epot = epot + sw*vterm
1020 +            
1021 +             a1 = 5.0d0 * ct_i * ct_j - ct_ij
1022 +            
1023 +             dudx = dudx + sw*pref*3.0d0*ri4 &
1024 +                             * (a1*xhat-ct_i*uz_j(1)-ct_j*uz_i(1))
1025 +             dudy = dudy + sw*pref*3.0d0*ri4 &
1026 +                             * (a1*yhat-ct_i*uz_j(2)-ct_j*uz_i(2))
1027 +             dudz = dudz + sw*pref*3.0d0*ri4 &
1028 +                             * (a1*zhat-ct_i*uz_j(3)-ct_j*uz_i(3))
1029 +            
1030 +             duduz_i(1) = duduz_i(1) + sw*pref*(ri3*(uz_j(1)-3.0d0*ct_j*xhat) &
1031 +                  - preRF2*uz_j(1))
1032 +             duduz_i(2) = duduz_i(2) + sw*pref*(ri3*(uz_j(2)-3.0d0*ct_j*yhat) &
1033 +                  - preRF2*uz_j(2))
1034 +             duduz_i(3) = duduz_i(3) + sw*pref*(ri3*(uz_j(3)-3.0d0*ct_j*zhat) &
1035 +                  - preRF2*uz_j(3))
1036 +             duduz_j(1) = duduz_j(1) + sw*pref*(ri3*(uz_i(1)-3.0d0*ct_i*xhat) &
1037 +                  - preRF2*uz_i(1))
1038 +             duduz_j(2) = duduz_j(2) + sw*pref*(ri3*(uz_i(2)-3.0d0*ct_i*yhat) &
1039 +                  - preRF2*uz_i(2))
1040 +             duduz_j(3) = duduz_j(3) + sw*pref*(ri3*(uz_i(3)-3.0d0*ct_i*zhat) &
1041 +                  - preRF2*uz_i(3))
1042  
1043            else
1044               if (i_is_SplitDipole) then
# Line 1080 | Line 1201 | contains
1201  
1202      if (do_pot) then
1203   #ifdef IS_MPI
1204 <       pot_row(atom1) = pot_row(atom1) + 0.5d0*epot
1205 <       pot_col(atom2) = pot_col(atom2) + 0.5d0*epot
1204 >       pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot
1205 >       pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot
1206   #else
1207         pot = pot + epot
1208   #endif
# Line 1186 | Line 1307 | contains
1307      return
1308    end subroutine doElectrostaticPair
1309  
1189 –  !! calculates the switching functions and their derivatives for a given
1190 –  subroutine calc_switch(r, mu, scale, dscale)
1191 –
1192 –    real (kind=dp), intent(in) :: r, mu
1193 –    real (kind=dp), intent(inout) :: scale, dscale
1194 –    real (kind=dp) :: rl, ru, mulow, minRatio, temp, scaleVal
1195 –
1196 –    ! distances must be in angstroms
1197 –    rl = 2.75d0
1198 –    ru = 3.75d0
1199 –    mulow = 0.0d0 !3.3856d0 ! 1.84 * 1.84
1200 –    minRatio = mulow / (mu*mu)
1201 –    scaleVal = 1.0d0 - minRatio
1202 –    
1203 –    if (r.lt.rl) then
1204 –       scale = minRatio
1205 –       dscale = 0.0d0
1206 –    elseif (r.gt.ru) then
1207 –       scale = 1.0d0
1208 –       dscale = 0.0d0
1209 –    else
1210 –       scale = 1.0d0 - scaleVal*((ru + 2.0d0*r - 3.0d0*rl) * (ru-r)**2) &
1211 –                        / ((ru - rl)**3)
1212 –       dscale = -scaleVal * 6.0d0 * (r-ru)*(r-rl)/((ru - rl)**3)    
1213 –    endif
1214 –        
1215 –    return
1216 –  end subroutine calc_switch
1217 –
1310    subroutine destroyElectrostaticTypes()
1311  
1312      if(allocated(ElectrostaticMap)) deallocate(ElectrostaticMap)
1313  
1314    end subroutine destroyElectrostaticTypes
1315  
1316 +  subroutine rf_self_self(atom1, eFrame, rfpot, t, do_pot)
1317 +    logical, intent(in) :: do_pot
1318 +    integer, intent(in) :: atom1
1319 +    integer :: atid1
1320 +    real(kind=dp), dimension(9,nLocal) :: eFrame
1321 +    real(kind=dp), dimension(3,nLocal) :: t
1322 +    real(kind=dp) :: mu1
1323 +    real(kind=dp) :: preVal, epot, rfpot
1324 +    real(kind=dp) :: eix, eiy, eiz
1325 +
1326 +    ! this is a local only array, so we use the local atom type id's:
1327 +    atid1 = atid(atom1)
1328 +    
1329 +    if (ElectrostaticMap(atid1)%is_Dipole) then
1330 +       mu1 = getDipoleMoment(atid1)
1331 +      
1332 +       preVal = pre22 * preRF2 * mu1*mu1
1333 +       rfpot = rfpot - 0.5d0*preVal
1334 +
1335 +       ! The self-correction term adds into the reaction field vector
1336 +      
1337 +       eix = preVal * eFrame(3,atom1)
1338 +       eiy = preVal * eFrame(6,atom1)
1339 +       eiz = preVal * eFrame(9,atom1)
1340 +
1341 +       ! once again, this is self-self, so only the local arrays are needed
1342 +       ! even for MPI jobs:
1343 +      
1344 +       t(1,atom1)=t(1,atom1) - eFrame(6,atom1)*eiz + &
1345 +            eFrame(9,atom1)*eiy
1346 +       t(2,atom1)=t(2,atom1) - eFrame(9,atom1)*eix + &
1347 +            eFrame(3,atom1)*eiz
1348 +       t(3,atom1)=t(3,atom1) - eFrame(3,atom1)*eiy + &
1349 +            eFrame(6,atom1)*eix
1350 +
1351 +    endif
1352 +    
1353 +    return
1354 +  end subroutine rf_self_self
1355 +
1356   end module electrostatic_module

Diff Legend

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