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 |