| 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 |
