1 |
subroutine reaction_field(pot, pre, natoms, mu, ulx, uly, ulz, tlx, tly, tlz, & |
2 |
rflx, rfly, rflz, is_dipole ) |
3 |
|
4 |
!passed variables |
5 |
|
6 |
integer :: natoms ! the number of atoms |
7 |
|
8 |
real(kind=8) :: pot ! hash, reefer, etc. aka. potential energy |
9 |
real(kind=8) :: pre ! a pretem multiplier for the reaction field |
10 |
! pre converts from mu in units of debye to kcal/mol |
11 |
|
12 |
|
13 |
|
14 |
!passed arrays |
15 |
|
16 |
logical(kind=2), dimension(natoms) :: is_dipole ! dipole identities |
17 |
|
18 |
real(kind=8), dimension(natoms) :: mu ! the dipole moment |
19 |
real(kind=8), dimension(natoms) :: ulx, uly, ulz ! the unitvectors |
20 |
real(kind=8), dimension(natoms) :: tlx, tly, tlz ! the torques |
21 |
real(kind=8), dimension(natoms) :: rflx, rfly, rflz ! reaction field |
22 |
|
23 |
! local variables |
24 |
|
25 |
integer :: i ! counter variable |
26 |
|
27 |
real(kind=8) :: dot ! the dot product |
28 |
real(kind=8) :: dVdu1x, dVdu1y, dVdu1z ! derivative of Potential wrt u |
29 |
|
30 |
|
31 |
! do single loop to compute torques on dipoles: |
32 |
|
33 |
do i = 1, natoms |
34 |
if ( is_dipole(i) ) then |
35 |
|
36 |
! remember that the sum over the cavity is supposed to include |
37 |
! the dipole at the center: |
38 |
|
39 |
dVdu1x = rflx(i) + (0.5_8 * ulx(i) * mu(i)) |
40 |
dVdu1y = rfly(i) + (0.5_8 * uly(i) * mu(i)) |
41 |
dVdu1z = rflz(i) + (0.5_8 * ulz(i) * mu(i)) |
42 |
|
43 |
rflx(i) = rflx(i) + (ulx(i) * mu(i)) |
44 |
rfly(i) = rfly(i) + (uly(i) * mu(i)) |
45 |
rflz(i) = rflz(i) + (ulz(i) * mu(i)) |
46 |
|
47 |
! The torque contribution is dipole cross reaction_field |
48 |
|
49 |
tlx(i) = tlx(i) - & |
50 |
( pre * mu(i) * ( (uly(i) * dVdu1z) - (ulz(i) * dVdu1y) ) ) |
51 |
|
52 |
tly(i) = tly(i) - & |
53 |
( pre * mu(i) * ( (ulz(i) * dVdu1x) - (ulx(i) * dVdu1z) ) ) |
54 |
|
55 |
tlz(i) = tlz(i) - & |
56 |
( pre * mu(i) * ( (ulx(i) * dVdu1y) - (uly(i) * dVdu1x) ) ) |
57 |
|
58 |
! the potential contribution is - 0.5 * dipole dot reaction_field |
59 |
|
60 |
dot = (ulx(i) * rflx(i)) + (uly(i) * rfly(i)) + (ulz(i) * rflz(i)) |
61 |
|
62 |
pot = pot + (0.5_8 * pre * mu(i) * dot) |
63 |
|
64 |
endif |
65 |
|
66 |
enddo |
67 |
|
68 |
end subroutine reaction_field |