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 2722 by gezelter, Thu Apr 20 18:24:24 2006 UTC vs.
Revision 2724 by chrisfen, Fri Apr 21 03:19:52 2006 UTC

# Line 77 | Line 77 | module electrostatic_module
77    real(kind=dp), parameter :: pre14 = 69.13373d0
78  
79    real(kind=dp), parameter :: zero = 0.0d0
80 +  
81 +  !! number of points for electrostatic splines
82 +  integer, parameter :: np = 100
83  
84    !! variables to handle different summation methods for long-range
85    !! electrostatics:
# Line 113 | Line 116 | module electrostatic_module
116    real(kind=dp), save :: f2c = 0.0_DP
117    real(kind=dp), save :: f3c = 0.0_DP
118    real(kind=dp), save :: f4c = 0.0_DP
119 +  real(kind=dp), save :: df0 = 0.0_DP
120 +  type(cubicSpline), save :: f0spline
121 +  logical, save :: haveElectroSpline = .false.
122  
123 +
124   #if defined(__IFC) || defined(__PGI)
125   ! error function for ifc version > 7.
126    double precision, external :: derfc
# Line 124 | Line 131 | module electrostatic_module
131    public :: setElectrostaticCutoffRadius
132    public :: setDampingAlpha
133    public :: setReactionFieldDielectric
134 <  public :: buildElectroSplines
134 >  public :: buildElectroSpline
135    public :: newElectrostaticType
136    public :: setCharge
137    public :: setDipoleMoment
# Line 196 | Line 203 | contains
203      haveDielectric = .true.
204    end subroutine setReactionFieldDielectric
205  
206 <  subroutine buildElectroSplines()
207 <  end subroutine buildElectroSplines
206 >  subroutine buildElectroSpline()
207 >    real( kind = dp ), dimension(np) :: xvals, yvals
208 >    real( kind = dp ) :: dx, rmin, rval
209 >    integer :: i
210 >
211 >    rmin = 0.0d0
212 >
213 >    dx = (defaultCutoff-rmin) / dble(np-1)
214 >    
215 >    do i = 1, np
216 >       rval = rmin + dble(i-1)*dx
217 >       xvals(i) = rval
218 >       yvals(i) = derfc(dampingAlpha*rval)
219 >    enddo
220  
221 +    call newSpline(f0spline, xvals, yvals, .true.)
222 +
223 +    haveElectroSpline = .true.
224 +  end subroutine buildElectroSpline
225 +
226    subroutine newElectrostaticType(c_ident, is_Charge, is_Dipole, &
227         is_SplitDipole, is_Quadrupole, is_Tap, status)
228  
# Line 450 | Line 474 | contains
474        
475      endif
476  
477 +    if (.not.haveElectroSpline) then
478 +       call buildElectroSpline()
479 +    end if
480 +
481      summationMethodChecked = .true.
482    end subroutine checkSummationMethod
483  
# Line 653 | Line 681 | contains
681  
682         if (j_is_Charge) then
683            if (screeningMethod .eq. DAMPED) then
684 <             f0 = derfc(dampingAlpha*rij)
685 <             varEXP = exp(-alpha2*rij*rij)
686 <             f1 = alphaPi*rij*varEXP + f0
684 >             call lookupUniformSpline1d(f0spline, rij, f0, df0)
685 >             f1 = -rij * df0 + f0
686 > !!$             f0 = derfc(dampingAlpha*rij)
687 > !!$             varEXP = exp(-alpha2*rij*rij)
688 > !!$             f1 = alphaPi*rij*varEXP + f0
689            endif
690  
691            preVal = pre11 * q_i * q_j
# Line 695 | Line 725 | contains
725  
726         if (j_is_Dipole) then
727            if (screeningMethod .eq. DAMPED) then
728 <             f0 = derfc(dampingAlpha*rij)
729 <             varEXP = exp(-alpha2*rij*rij)
730 <             f1 = alphaPi*rij*varEXP + f0
731 <             f3 = alphaPi*2.0d0*alpha2*varEXP*rij*rij*rij
728 >             call lookupUniformSpline1d(f0spline, rij, f0, df0)
729 >             f1 = -rij * df0 + f0
730 >             f3 = -2.0d0*alpha2*df0*rij*rij*rij
731 > !!$             f0 = derfc(dampingAlpha*rij)
732 > !!$             varEXP = exp(-alpha2*rij*rij)
733 > !!$             f1 = alphaPi*rij*varEXP + f0
734 > !!$             f3 = alphaPi*2.0d0*alpha2*varEXP*rij*rij*rij
735            endif
736  
737            pref = pre12 * q_i * mu_j
# Line 764 | Line 797 | contains
797  
798         if (j_is_Quadrupole) then
799            if (screeningMethod .eq. DAMPED) then
800 <             f0 = derfc(dampingAlpha*rij)
801 <             varEXP = exp(-alpha2*rij*rij)
802 <             f1 = alphaPi*rij*varEXP + f0
803 <             f2 = alphaPi*2.0d0*alpha2*varEXP
800 >             call lookupUniformSpline1d(f0spline, rij, f0, df0)
801 > !!$             f0 = derfc(dampingAlpha*rij)
802 > !!$             varEXP = exp(-alpha2*rij*rij)
803 > !!$             f1 = alphaPi*rij*varEXP + f0
804 > !!$             f2 = alphaPi*2.0d0*alpha2*varEXP
805 >             f1 = -rij * df0 + f0
806 >             f2 = -2.0d0*alpha2*df0
807               f3 = f2*rij*rij*rij
808               f4 = 2.0d0*alpha2*f2*rij
809            endif
# Line 834 | Line 870 | contains
870  
871         if (j_is_Charge) then
872            if (screeningMethod .eq. DAMPED) then
873 <             f0 = derfc(dampingAlpha*rij)
874 <             varEXP = exp(-alpha2*rij*rij)
875 <             f1 = alphaPi*rij*varEXP + f0
876 <             f3 = alphaPi*2.0d0*alpha2*varEXP*rij*rij*rij
873 >             call lookupUniformSpline1d(f0spline, rij, f0, df0)
874 >             f1 = -rij * df0 + f0
875 >             f3 = -2.0d0*alpha2*df0*rij*rij*rij
876 > !!$             f0 = derfc(dampingAlpha*rij)
877 > !!$             varEXP = exp(-alpha2*rij*rij)
878 > !!$             f1 = alphaPi*rij*varEXP + f0
879 > !!$             f3 = alphaPi*2.0d0*alpha2*varEXP*rij*rij*rij
880            endif
881            
882            pref = pre12 * q_j * mu_i
# Line 934 | Line 973 | contains
973        
974         if (j_is_Dipole) then
975            if (screeningMethod .eq. DAMPED) then
976 <             f0 = derfc(dampingAlpha*rij)
977 <             varEXP = exp(-alpha2*rij*rij)
978 <             f1 = alphaPi*rij*varEXP + f0
979 <             f2 = alphaPi*2.0d0*alpha2*varEXP
976 >             call lookupUniformSpline1d(f0spline, rij, f0, df0)
977 > !!$             f0 = derfc(dampingAlpha*rij)
978 > !!$             varEXP = exp(-alpha2*rij*rij)
979 > !!$             f1 = alphaPi*rij*varEXP + f0
980 > !!$             f2 = alphaPi*2.0d0*alpha2*varEXP
981 >             f1 = -rij * df0 + f0
982 >             f2 = -2.0d0*alpha2*df0
983               f3 = f2*rij*rij*rij
984               f4 = 2.0d0*alpha2*f3*rij*rij
985            endif
# Line 1054 | Line 1096 | contains
1096      if (i_is_Quadrupole) then
1097         if (j_is_Charge) then
1098            if (screeningMethod .eq. DAMPED) then
1099 <             f0 = derfc(dampingAlpha*rij)
1100 <             varEXP = exp(-alpha2*rij*rij)
1101 <             f1 = alphaPi*rij*varEXP + f0
1102 <             f2 = alphaPi*2.0d0*alpha2*varEXP
1099 >             call lookupUniformSpline1d(f0spline, rij, f0, df0)
1100 > !!$             f0 = derfc(dampingAlpha*rij)
1101 > !!$             varEXP = exp(-alpha2*rij*rij)
1102 > !!$             f1 = alphaPi*rij*varEXP + f0
1103 > !!$             f2 = alphaPi*2.0d0*alpha2*varEXP
1104 >             f1 = -rij * df0 + f0
1105 >             f2 = -2.0d0*alpha2*df0
1106               f3 = f2*rij*rij*rij
1107               f4 = 2.0d0*alpha2*f2*rij
1108            endif

Diff Legend

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