| 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 |
< |
|
| 93 |
< |
|
| 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 |
| 96 |
> |
#endif |
| 97 |
> |
|
| 98 |
|
public :: setElectrostaticSummationMethod |
| 99 |
|
public :: setElectrostaticCutoffRadius |
| 100 |
|
public :: setDampedWolfAlpha |
| 136 |
|
|
| 137 |
|
summationMethod = the_ESM |
| 138 |
|
|
| 135 |
– |
if (summationMethod .eq. UNDAMPED_WOLF) is_Undamped_Wolf = .true. |
| 136 |
– |
if (summationMethod .eq. DAMPED_WOLF) is_Damped_Wolf = .true. |
| 139 |
|
end subroutine setElectrostaticSummationMethod |
| 140 |
|
|
| 141 |
|
subroutine setElectrostaticCutoffRadius(thisRcut) |
| 389 |
|
endif |
| 390 |
|
|
| 391 |
|
constEXP = exp(-dampingAlpha*dampingAlpha*defaultCutoff*defaultCutoff) |
| 392 |
< |
constERFC = erfc(dampingAlpha*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) |
| 641 |
+ |
dudz = dudz + dudr * d(3) |
| 642 |
+ |
|
| 643 |
|
else |
| 644 |
|
|
| 645 |
|
vterm = pre11 * q_i * q_j * riji |
