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: |
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 |
131 |
|
public :: setElectrostaticCutoffRadius |
132 |
|
public :: setDampingAlpha |
133 |
|
public :: setReactionFieldDielectric |
134 |
< |
public :: buildElectroSplines |
134 |
> |
public :: buildElectroSpline |
135 |
|
public :: newElectrostaticType |
136 |
|
public :: setCharge |
137 |
|
public :: setDipoleMoment |
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 |
|
|
474 |
|
|
475 |
|
endif |
476 |
|
|
477 |
+ |
if (.not.haveElectroSpline) then |
478 |
+ |
call buildElectroSpline() |
479 |
+ |
end if |
480 |
+ |
|
481 |
|
summationMethodChecked = .true. |
482 |
|
end subroutine checkSummationMethod |
483 |
|
|
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 |
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 |
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 |
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 |
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 |
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 |