| 1 |
subroutine dipole_rf( atom1, atom2, rrfsq, r2, pot, rxij, ryij, rzij, rrf, & |
| 2 |
dielect, rt, natoms, ulx, uly, ulz, mu, flx, fly, flz, tlx, tly, tlz, & |
| 3 |
rflx, rfly, rflz, prerf ) |
| 4 |
|
| 5 |
!passed arguments |
| 6 |
|
| 7 |
integer :: atom1, atom2 ! atom i and j |
| 8 |
integer :: natoms ! the number of atoms in the system |
| 9 |
|
| 10 |
double precision rrfsq ! the square of the reaction field radius |
| 11 |
double precision r2 ! the square of the distance twixt i and j |
| 12 |
double precision pot ! hash, reefer, aka. the potential energy |
| 13 |
double precision rxij, ryij, rzij ! the rij vector components |
| 14 |
double precision rrf ! the reaction field radius |
| 15 |
double precision dielect ! the dielectric constant of the reaction field |
| 16 |
double precision rt ! the radius of the reaction field taper |
| 17 |
double precision prerf ! dipole equation prefactor |
| 18 |
|
| 19 |
!passed arrays |
| 20 |
|
| 21 |
double precision, dimension(natoms) :: ulx, uly, ulz ! lab frame unit vector |
| 22 |
double precision, dimension(natoms) :: mu ! the dipole moments |
| 23 |
double precision, dimension(natoms) :: flx, fly, flz ! forces |
| 24 |
double precision, dimension(natoms) :: tlx, tly, tlz ! torques |
| 25 |
double precision, dimension(natoms) :: rflx, rfly, rflz ! reaction field |
| 26 |
|
| 27 |
!local variables |
| 28 |
|
| 29 |
double precision dx, dy, dz ! the rji components |
| 30 |
double precision rij ! the distance twixt i and j |
| 31 |
double precision dfact1, dfact2, dip2, pre ! assorted prefactors |
| 32 |
double precision r3, r5 ! the cube and 5th of rij |
| 33 |
double precision dudx, dudy, dudz ! derivative of u wrt x, y, and z |
| 34 |
double precision dudu1x, dudu1y, dudu1z ! derivative of u wrt u of atom i |
| 35 |
double precision dudu2x, dudu2y, dudu2z ! derivative of u wrt u oj atom j |
| 36 |
double precision rdotu1, rdotu2 ! r dot u of i and j |
| 37 |
double precision u1dotu2 ! u of i dot u of j |
| 38 |
double precision taper, dtdr ! taper variables |
| 39 |
double precision vterm ! dipole potential term |
| 40 |
|
| 41 |
double precision taper_2, dtdr_2 |
| 42 |
logical mine |
| 43 |
|
| 44 |
dx = rxij |
| 45 |
dy = ryij |
| 46 |
dz = rzij |
| 47 |
|
| 48 |
mine = .true. |
| 49 |
|
| 50 |
! pre converts from mu in units of debye to kcal/mol |
| 51 |
|
| 52 |
pre = 14.38362d0 |
| 53 |
rij = dsqrt( r2 ) |
| 54 |
|
| 55 |
if ( rij .le. rrf ) then |
| 56 |
|
| 57 |
! cubic taper |
| 58 |
if ( rij .lt. rt ) then |
| 59 |
taper = 1.0d0 |
| 60 |
dtdr = 0.0d0 |
| 61 |
else |
| 62 |
taper = ( rrf + (2.0d0 * rij) - (3.0d0 * rt) ) * ( rrf - rij )**2 / & |
| 63 |
( ( rrf - rt )**3 ) |
| 64 |
dtdr = 6.0d0 * ( (rij * rij) - (rij * rt) - (rij * rrf) + & |
| 65 |
(rrf * rt) ) / ( (rrf - rt)**3 ) |
| 66 |
|
| 67 |
endif |
| 68 |
|
| 69 |
r3 = r2 * rij |
| 70 |
r5 = r3 * r2 |
| 71 |
|
| 72 |
rdotu1 = (dx * ulx(atom1)) + (dy * uly(atom1)) + (dz * ulz(atom1)) |
| 73 |
rdotu2 = (dx * ulx(atom2)) + (dy * uly(atom2)) + (dz * ulz(atom2)) |
| 74 |
u1dotu2 = (ulx(atom1) * ulx(atom2)) + (uly(atom1) * uly(atom2)) + & |
| 75 |
(ulz(atom1) * ulz(atom2)) |
| 76 |
|
| 77 |
dip2 = pre * mu(atom1) * mu(atom2) |
| 78 |
|
| 79 |
dfact1 = 3.0d0 * dip2 / r2 |
| 80 |
dfact2 = 3.0d0 * dip2 / r5 |
| 81 |
|
| 82 |
vterm = dip2 * ( (u1dotu2 / r3) - 3.0d0 * (rdotu1 * rdotu2 / r5) ) |
| 83 |
|
| 84 |
pot = pot + vterm * taper |
| 85 |
|
| 86 |
dudx = ( -dfact1 * dx * & |
| 87 |
( (u1dotu2 / r3) - ( 5.0d0 * (rdotu1 * rdotu2) / r5 ) ) - & |
| 88 |
dfact2 * ( (ulx(atom1) * rdotu2) + (ulx(atom2) * rdotu1) ) )* & |
| 89 |
taper + & |
| 90 |
vterm * dtdr * dx / rij |
| 91 |
|
| 92 |
dudy = ( -dfact1 * dy * & |
| 93 |
( (u1dotu2 / r3) - ( 5.0d0 * (rdotu1 * rdotu2) / r5 ) ) - & |
| 94 |
dfact2 * ( (uly(atom1) * rdotu2) + (uly(atom2) * rdotu1) ) ) * & |
| 95 |
taper + & |
| 96 |
vterm * dtdr * dy / rij |
| 97 |
|
| 98 |
dudz = ( -dfact1 * dz * & |
| 99 |
( (u1dotu2 / r3) - ( 5.0d0 * (rdotu1 * rdotu2) / r5 ) ) - & |
| 100 |
dfact2 * ( (ulz(atom1) * rdotu2) + (ulz(atom2) * rdotu1) ) ) * & |
| 101 |
taper + & |
| 102 |
vterm * dtdr * dz / rij |
| 103 |
|
| 104 |
dudu1x = ( dip2 * ( (ulx(atom2) / r3) - ( 3.0d0 * dx * rdotu2 / r5 ) ) ) & |
| 105 |
* taper |
| 106 |
dudu1y = ( dip2 * ( (uly(atom2) / r3) - ( 3.0d0 * dy * rdotu2 / r5 ) ) ) & |
| 107 |
* taper |
| 108 |
dudu1z = ( dip2 * ( (ulz(atom2) / r3) - ( 3.0d0 * dz * rdotu2 / r5 ) ) ) & |
| 109 |
* taper |
| 110 |
|
| 111 |
dudu2x = ( dip2 * ( (ulx(atom1) / r3) - ( 3.0d0 * dx * rdotu1 / r5 ) ) ) & |
| 112 |
* taper |
| 113 |
dudu2y = ( dip2 * ( (uly(atom1) / r3) - ( 3.0d0 * dy * rdotu1 / r5 ) ) ) & |
| 114 |
* taper |
| 115 |
dudu2z = ( dip2 * ( (ulz(atom1) / r3) - ( 3.0d0 * dz * rdotu1 / r5 ) ) ) & |
| 116 |
* taper |
| 117 |
|
| 118 |
flx(atom1) = flx(atom1) + dudx |
| 119 |
fly(atom1) = fly(atom1) + dudy |
| 120 |
flz(atom1) = flz(atom1) + dudz |
| 121 |
|
| 122 |
flx(atom2) = flx(atom2) - dudx |
| 123 |
fly(atom2) = fly(atom2) - dudy |
| 124 |
flz(atom2) = flz(atom2) - dudz |
| 125 |
|
| 126 |
tlx(atom1) = tlx(atom1) - (uly(atom1) * dudu1z) + (ulz(atom1) * dudu1y) |
| 127 |
tly(atom1) = tly(atom1) - (ulz(atom1) * dudu1x) + (ulx(atom1) * dudu1z) |
| 128 |
tlz(atom1) = tlz(atom1) - (ulx(atom1) * dudu1y) + (uly(atom1) * dudu1x) |
| 129 |
|
| 130 |
tlx(atom2) = tlx(atom2) - (uly(atom2) * dudu2z) + (ulz(atom2) * dudu2y) |
| 131 |
tly(atom2) = tly(atom2) - (ulz(atom2) * dudu2x) + (ulx(atom2) * dudu2z) |
| 132 |
tlz(atom2) = tlz(atom2) - (ulx(atom2) * dudu2y) + (uly(atom2) * dudu2x) |
| 133 |
|
| 134 |
endif |
| 135 |
|
| 136 |
return |
| 137 |
end subroutine dipole_rf |
