| 1 |
subroutine dipole_dipole(atom1, atom2, dx, dy, dz, rij, pot) |
| 2 |
|
| 3 |
include 'sizes.inc' |
| 4 |
include 'simulation.inc' |
| 5 |
|
| 6 |
integer atom1, atom2 |
| 7 |
double precision dx, dy, dz, rij |
| 8 |
double precision dfact1, dfact2, dip2, r2, r3, r5, pre |
| 9 |
double precision dudx, dudy, dudz, dudu1x, dudu1y, dudu1z |
| 10 |
double precision dudu2x, dudu2y, dudu2z, rdotu1, rdotu2, u1dotu2 |
| 11 |
double precision taper, dtdr, vterm |
| 12 |
|
| 13 |
! pre converts from mu in units of debye to kcal/mol |
| 14 |
|
| 15 |
pre = 14.38362d0 |
| 16 |
|
| 17 |
if (rij.le.rrf) then |
| 18 |
|
| 19 |
if (rij.lt.rt) then |
| 20 |
taper = 1.0d0 |
| 21 |
dtdr = 0.0d0 |
| 22 |
else |
| 23 |
taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3) |
| 24 |
dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3) |
| 25 |
endif |
| 26 |
|
| 27 |
r2 = rij*rij |
| 28 |
r3 = r2*rij |
| 29 |
r5 = r3*r2 |
| 30 |
|
| 31 |
rdotu1 = dx*ulx(atom1) + dy*uly(atom1) + dz*ulz(atom1) |
| 32 |
rdotu2 = dx*ulx(atom2) + dy*uly(atom2) + dz*ulz(atom2) |
| 33 |
u1dotu2 = ulx(atom1)*ulx(atom2) + uly(atom1)*uly(atom2) + & |
| 34 |
ulz(atom1)*ulz(atom2) |
| 35 |
|
| 36 |
dip2 = pre*mu(atom1)*mu(atom2) |
| 37 |
|
| 38 |
dfact1 = 3.0d0*dip2 / r2 |
| 39 |
dfact2 = 3.0d0*dip2 / r5 |
| 40 |
|
| 41 |
vterm = dip2*((u1dotu2/r3) - 3.0d0*(rdotu1*rdotu2/r5)) |
| 42 |
|
| 43 |
pot = pot + vterm*taper |
| 44 |
|
| 45 |
dudx = (-dfact1 * dx * ((u1dotu2/r3) - & |
| 46 |
(5.0d0*(rdotu1*rdotu2)/r5)) - & |
| 47 |
dfact2*(ulx(atom1)*rdotu2 + ulx(atom2)*rdotu1))*taper + & |
| 48 |
vterm*dtdr*dx/rij |
| 49 |
|
| 50 |
dudy = (-dfact1 * dy * ((u1dotu2/r3) - & |
| 51 |
(5.0d0*(rdotu1*rdotu2)/r5)) - & |
| 52 |
dfact2*(uly(atom1)*rdotu2 + uly(atom2)*rdotu1))*taper + & |
| 53 |
vterm*dtdr*dy/rij |
| 54 |
|
| 55 |
dudz = (-dfact1 * dz * ((u1dotu2/r3) - & |
| 56 |
(5.0d0*(rdotu1*rdotu2)/r5)) - & |
| 57 |
dfact2*(ulz(atom1)*rdotu2 + ulz(atom2)*rdotu1))*taper + & |
| 58 |
vterm*dtdr*dz/rij |
| 59 |
|
| 60 |
dudu1x = (dip2*((ulx(atom2)/r3) - (3.0d0*dx*rdotu2/r5)))*taper |
| 61 |
dudu1y = (dip2*((uly(atom2)/r3) - (3.0d0*dy*rdotu2/r5)))*taper |
| 62 |
dudu1z = (dip2*((ulz(atom2)/r3) - (3.0d0*dz*rdotu2/r5)))*taper |
| 63 |
|
| 64 |
dudu2x = (dip2*((ulx(atom1)/r3) - (3.0d0*dx*rdotu1/r5)))*taper |
| 65 |
dudu2y = (dip2*((uly(atom1)/r3) - (3.0d0*dy*rdotu1/r5)))*taper |
| 66 |
dudu2z = (dip2*((ulz(atom1)/r3) - (3.0d0*dz*rdotu1/r5)))*taper |
| 67 |
|
| 68 |
flx(atom1) = flx(atom1) + dudx |
| 69 |
fly(atom1) = fly(atom1) + dudy |
| 70 |
flz(atom1) = flz(atom1) + dudz |
| 71 |
|
| 72 |
flx(atom2) = flx(atom2) - dudx |
| 73 |
fly(atom2) = fly(atom2) - dudy |
| 74 |
flz(atom2) = flz(atom2) - dudz |
| 75 |
|
| 76 |
tlx(atom1) = tlx(atom1) - uly(atom1)*dudu1z + ulz(atom1)*dudu1y |
| 77 |
tly(atom1) = tly(atom1) - ulz(atom1)*dudu1x + ulx(atom1)*dudu1z |
| 78 |
tlz(atom1) = tlz(atom1) - ulx(atom1)*dudu1y + uly(atom1)*dudu1x |
| 79 |
|
| 80 |
tlx(atom2) = tlx(atom2) - uly(atom2)*dudu2z + ulz(atom2)*dudu2y |
| 81 |
tly(atom2) = tly(atom2) - ulz(atom2)*dudu2x + ulx(atom2)*dudu2z |
| 82 |
tlz(atom2) = tlz(atom2) - ulx(atom2)*dudu2y + uly(atom2)*dudu2x |
| 83 |
|
| 84 |
virial = virial + ( dx*dudx + dy*dudy + dz*dudz ) |
| 85 |
|
| 86 |
endif |
| 87 |
|
| 88 |
return |
| 89 |
end subroutine dipole_dipole |
