87 |
|
real(kind=dp), save :: rcuti2 = 0.0_dp |
88 |
|
real(kind=dp), save :: rcuti3 = 0.0_dp |
89 |
|
real(kind=dp), save :: rcuti4 = 0.0_dp |
90 |
< |
logical, save :: is_Undamped_Wolf = .false. |
91 |
< |
logical, save :: is_Damped_Wolf = .false. |
92 |
< |
|
90 |
> |
real(kind=dp), save :: alphaPi = 0.0_dp |
91 |
> |
real(kind=dp), save :: invRootPi = 0.0_dp |
92 |
> |
|
93 |
|
#ifdef __IFC |
94 |
|
! error function for ifc version > 7. |
95 |
|
double precision, external :: derfc |
136 |
|
|
137 |
|
summationMethod = the_ESM |
138 |
|
|
139 |
– |
if (summationMethod .eq. UNDAMPED_WOLF) is_Undamped_Wolf = .true. |
140 |
– |
if (summationMethod .eq. DAMPED_WOLF) is_Damped_Wolf = .true. |
139 |
|
end subroutine setElectrostaticSummationMethod |
140 |
|
|
141 |
|
subroutine setElectrostaticCutoffRadius(thisRcut) |
390 |
|
|
391 |
|
constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff) |
392 |
|
constERFC = derfc(dampingAlpha*defaultCutoff) |
393 |
+ |
invRootPi = 0.56418958354775628695d0 |
394 |
+ |
alphaPi = 2*dampingAlpha*invRootPi |
395 |
|
|
396 |
|
haveDWAconstants = .true. |
397 |
|
endif |
447 |
|
real (kind=dp) :: xhat, yhat, zhat |
448 |
|
real (kind=dp) :: dudx, dudy, dudz |
449 |
|
real (kind=dp) :: scale, sc2, bigR |
450 |
+ |
real (kind=dp) :: varERFC, varEXP |
451 |
|
|
452 |
|
if (.not.allocated(ElectrostaticMap)) then |
453 |
|
call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!") |
617 |
|
vpair = vpair + vterm |
618 |
|
epot = epot + sw*vterm |
619 |
|
|
620 |
< |
dudr = - sw * pre11 * q_i * q_j * (riji*riji*riji - rcuti2*rcuti) |
620 |
> |
dudr = -sw*pre11*q_i*q_j * (riji*riji*riji - rcuti2*rcuti) |
621 |
> |
|
622 |
> |
dudx = dudx + dudr * d(1) |
623 |
> |
dudy = dudy + dudr * d(2) |
624 |
> |
dudz = dudz + dudr * d(3) |
625 |
> |
|
626 |
> |
elseif (summationMethod .eq. DAMPED_WOLF) then |
627 |
> |
|
628 |
> |
varERFC = derfc(dampingAlpha*rij) |
629 |
> |
varEXP = exp(-dampingAlpha*dampingAlpha*rij*rij) |
630 |
> |
vterm = pre11 * q_i * q_j * (varERFC*riji - constERFC*rcuti) |
631 |
> |
vpair = vpair + vterm |
632 |
> |
epot = epot + sw*vterm |
633 |
> |
|
634 |
> |
dudr = -sw*pre11*q_i*q_j * ( riji*(varERFC*riji*riji & |
635 |
> |
+ alphaPi*varEXP) & |
636 |
> |
- rcuti*(constERFC*rcuti2 & |
637 |
> |
+ alphaPi*constEXP) ) |
638 |
|
|
639 |
|
dudx = dudx + dudr * d(1) |
640 |
|
dudy = dudy + dudr * d(2) |