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 2717 by gezelter, Mon Apr 17 21:49:12 2006 UTC vs.
Revision 2724 by chrisfen, Fri Apr 21 03:19:52 2006 UTC

# Line 76 | Line 76 | module electrostatic_module
76    !! This unit is also known affectionately as an esu centi-barn.
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:
86    integer, save :: summationMethod = NONE
# Line 111 | 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 122 | 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 194 | 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 448 | 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 550 | Line 580 | contains
580         if (i_is_SplitDipole) then
581            d_i = ElectrostaticMap(me1)%split_dipole_distance
582         endif
583 <
583 >       duduz_i = zero
584      endif
585  
586      if (i_is_Quadrupole) then
# Line 581 | Line 611 | contains
611         cx_i = ux_i(1)*xhat + ux_i(2)*yhat + ux_i(3)*zhat
612         cy_i = uy_i(1)*xhat + uy_i(2)*yhat + uy_i(3)*zhat
613         cz_i = uz_i(1)*xhat + uz_i(2)*yhat + uz_i(3)*zhat
614 +       dudux_i = zero
615 +       duduy_i = zero
616 +       duduz_i = zero
617      endif
618  
619      if (j_is_Charge) then
# Line 603 | Line 636 | contains
636         if (j_is_SplitDipole) then
637            d_j = ElectrostaticMap(me2)%split_dipole_distance
638         endif
639 +       duduz_j = zero
640      endif
641  
642      if (j_is_Quadrupole) then
# Line 633 | Line 667 | contains
667         cx_j = ux_j(1)*xhat + ux_j(2)*yhat + ux_j(3)*zhat
668         cy_j = uy_j(1)*xhat + uy_j(2)*yhat + uy_j(3)*zhat
669         cz_j = uz_j(1)*xhat + uz_j(2)*yhat + uz_j(3)*zhat
670 +       dudux_j = zero
671 +       duduy_j = zero
672 +       duduz_j = zero
673      endif
674    
675 <    epot = 0.0d0
676 <    dudx = 0.0d0
677 <    dudy = 0.0d0
678 <    dudz = 0.0d0
675 >    epot = zero
676 >    dudx = zero
677 >    dudy = zero
678 >    dudz = zero  
679  
643 –    dudux_i = 0.0d0
644 –    duduy_i = 0.0d0
645 –    duduz_i = 0.0d0
646 –
647 –    dudux_j = 0.0d0
648 –    duduy_j = 0.0d0
649 –    duduz_j = 0.0d0
650 –
680      if (i_is_Charge) then
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 694 | 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 763 | 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 833 | 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 933 | 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 1053 | 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
# Line 1326 | Line 1372 | contains
1372         call checkSummationMethod()
1373      endif
1374  
1375 <    dudx = 0.0d0
1376 <    dudy = 0.0d0
1377 <    dudz = 0.0d0
1375 >    dudx = zero
1376 >    dudy = zero
1377 >    dudz = zero
1378  
1379      riji = 1.0d0/rij
1380  

Diff Legend

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