1 |
mmeineke |
10 |
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 |