54 |
|
|
55 |
|
PRIVATE |
56 |
|
|
57 |
+ |
|
58 |
|
#define __FORTRAN90 |
59 |
+ |
#include "UseTheForce/DarkSide/fInteractionMap.h" |
60 |
|
#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" |
61 |
+ |
|
62 |
|
|
63 |
|
!! these prefactors convert the multipole interactions into kcal / mol |
64 |
|
!! all were computed assuming distances are measured in angstroms |
395 |
|
constERFC = derfc(dampingAlpha*defaultCutoff) |
396 |
|
invRootPi = 0.56418958354775628695d0 |
397 |
|
alphaPi = 2*dampingAlpha*invRootPi |
398 |
< |
|
398 |
> |
|
399 |
|
haveDWAconstants = .true. |
400 |
|
endif |
401 |
|
endif |
451 |
|
real (kind=dp) :: dudx, dudy, dudz |
452 |
|
real (kind=dp) :: scale, sc2, bigR |
453 |
|
real (kind=dp) :: varERFC, varEXP |
454 |
+ |
real (kind=dp) :: limScale |
455 |
|
|
456 |
|
if (.not.allocated(ElectrostaticMap)) then |
457 |
|
call handleError("electrostatic", "no ElectrostaticMap was present before first call of do_electrostatic_pair!") |
475 |
|
!! some variables we'll need independent of electrostatic type: |
476 |
|
|
477 |
|
riji = 1.0d0 / rij |
478 |
< |
|
478 |
> |
|
479 |
|
xhat = d(1) * riji |
480 |
|
yhat = d(2) * riji |
481 |
|
zhat = d(3) * riji |
621 |
|
vpair = vpair + vterm |
622 |
|
epot = epot + sw*vterm |
623 |
|
|
624 |
< |
dudr = -sw*pre11*q_i*q_j * (riji*riji*riji - rcuti2*rcuti) |
624 |
> |
dudr = -sw*pre11*q_i*q_j * (riji*riji-rcuti2)*riji |
625 |
|
|
626 |
|
dudx = dudx + dudr * d(1) |
627 |
|
dudy = dudy + dudr * d(2) |
635 |
|
vpair = vpair + vterm |
636 |
|
epot = epot + sw*vterm |
637 |
|
|
638 |
< |
dudr = -sw*pre11*q_i*q_j * ( riji*(varERFC*riji*riji & |
639 |
< |
+ alphaPi*varEXP) & |
640 |
< |
- rcuti*(constERFC*rcuti2 & |
641 |
< |
+ alphaPi*constEXP) ) |
638 |
> |
dudr = -sw*pre11*q_i*q_j * ( riji*((varERFC*riji*riji & |
639 |
> |
+ alphaPi*varEXP) & |
640 |
> |
- (constERFC*rcuti2 & |
641 |
> |
+ alphaPi*constEXP)) ) |
642 |
|
|
643 |
|
dudx = dudx + dudr * d(1) |
644 |
|
dudy = dudy + dudr * d(2) |
1084 |
|
|
1085 |
|
if (do_pot) then |
1086 |
|
#ifdef IS_MPI |
1087 |
< |
pot_row(atom1) = pot_row(atom1) + 0.5d0*epot |
1088 |
< |
pot_col(atom2) = pot_col(atom2) + 0.5d0*epot |
1087 |
> |
pot_row(ELECTROSTATIC_POT,atom1) = pot_row(ELECTROSTATIC_POT,atom1) + 0.5d0*epot |
1088 |
> |
pot_col(ELECTROSTATIC_POT,atom2) = pot_col(ELECTROSTATIC_POT,atom2) + 0.5d0*epot |
1089 |
|
#else |
1090 |
|
pot = pot + epot |
1091 |
|
#endif |