--- trunk/OOPSE_old/src/mdtools/libmdCode/calc_reaction_field.F90 2003/03/10 19:37:48 307 +++ trunk/OOPSE_old/src/mdtools/libmdCode/calc_reaction_field.F90 2003/03/18 19:09:12 361 @@ -1,141 +1,248 @@ -subroutine accumulate_rf(atom1, atom2, rij) +module reaction_field + use force_globals + use definitions + use atype_module + use vector_class +#ifdef IS_MPI + use mpiSimulation +#endif + implicit none - include 'sizes.inc' - include 'simulation.inc' + real(kind=dp), save :: rrf + real(kind=dp), save :: rt + real(kind=dp), save :: dielect + real(kind=dp), save :: rrfsq + real(kind=dp), save :: pre + logical, save :: rf_initialized = .false. - integer atom1, atom2 - double precision taper, rij - - if (rij.le.rrf) then - - if (rij.lt.rt) then - taper = 1.0d0 - else - taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3) - endif - - rflx(atom1) = rflx(atom1) + ulx(atom2)*mu(atom2)*taper - rfly(atom1) = rfly(atom1) + uly(atom2)*mu(atom2)*taper - rflz(atom1) = rflz(atom1) + ulz(atom2)*mu(atom2)*taper - - rflx(atom2) = rflx(atom2) + ulx(atom1)*mu(atom1)*taper - rfly(atom2) = rfly(atom2) + uly(atom1)*mu(atom1)*taper - rflz(atom2) = rflz(atom2) + ulz(atom1)*mu(atom1)*taper - - endif - return +contains -end subroutine accumulate_rf + subroutine initialize_rf(this_rrf, this_rt, this_dielect) + real(kind=dp), intent(in) :: this_rrf, this_rt, this_dielect -subroutine accumulate_self_rf() - - include 'sizes.inc' - include 'simulation.inc' - - integer i, ia, a1, atype1 - logical is_dipole_atype - external is_dipole_atype - - do i = 1, nmol - do ia = 1, na(i) - a1 = atom_index(i,ia) - - atype1 = atype(a1) - - if (is_dipole_atype(atype1)) then - - rflx(a1) = rflx(a1) + ulx(a1)*mu(a1) - rfly(a1) = rfly(a1) + uly(a1)*mu(a1) - rflz(a1) = rflz(a1) + ulz(a1)*mu(a1) - - endif - enddo - enddo - - return -end subroutine accumulate_self_rf - -subroutine reaction_field(pot) + rrf = this_rrf + rt = this_rt + dielect = this_dielect + + rrfsq = rrf * rrf + pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) + + rf_initialized = .true. + + return + end subroutine initialize_rf - include 'sizes.inc' - include 'simulation.inc' - - double precision rrfsq, pre - integer i, ia, a1, atype1 - logical is_dipole_atype - external is_dipole_atype + subroutine accumulate_rf(atom1, atom2, rij, u_l) - ! do single loop to compute torques on dipoles: - ! pre converts from mu in units of debye to kcal/mol + integer, intent(in) :: atom1, atom2 + real (kind = dp), intent(in) :: rij + real (kind = dp), dimension(:,:) :: u_l - rrfsq = rrf * rrf - pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) - - do i = 1, nmol - do ia = 1, na(i) - a1 = atom_index(i,ia) + integer :: me1, me2 + real (kind = dp) :: taper, mu1, mu2 + real (kind = dp), dimension(3) :: ul1 + real (kind = dp), dimension(3) :: ul2 - atype1 = atype(a1) - - if (is_dipole_atype(atype1)) then - - ! The torque contribution is dipole cross reaction_field - - tlx(a1) = tlx(a1) + pre*mu(a1)*(uly(a1)*rflz(a1) - ulz(a1)*rfly(a1)) - tly(a1) = tly(a1) + pre*mu(a1)*(ulz(a1)*rflx(a1) - ulx(a1)*rflz(a1)) - tlz(a1) = tlz(a1) + pre*mu(a1)*(ulx(a1)*rfly(a1) - uly(a1)*rflx(a1)) - - ! the potential contribution is -1/2 dipole dot reaction_field + if (.not.rf_initialized) then + write(default_error,*) 'Reaction field not initialized!' + return + endif + + if (rij.le.rrf) then - pot = pot - 0.5d0 * pre * mu(a1) * & - (rflx(a1)*ulx(a1) + rfly(a1)*uly(a1) + rflz(a1)*ulz(a1)) - - endif - - enddo - enddo + if (rij.lt.rt) then + taper = 1.0d0 + else + taper = (rrf + 2.0d0*rij - 3.0d0*rt)*(rrf-rij)**2/ ((rrf-rt)**3) + endif + +#ifdef IS_MPI + me1 = atid_Row(atom1) + ul1(1) = u_l_Row(1,atom1) + ul1(2) = u_l_Row(2,atom1) + ul1(3) = u_l_Row(3,atom1) + + me2 = atid_Col(atom2) + ul2(1) = u_l_Col(1,atom2) + ul2(2) = u_l_Col(2,atom2) + ul2(3) = u_l_Col(3,atom2) +#else + me1 = atid(atom1) + ul1(1) = u_l(1,atom1) + ul1(2) = u_l(2,atom1) + ul1(3) = u_l(3,atom1) + + me2 = atid(atom2) + ul2(1) = u_l(1,atom2) + ul2(2) = u_l(2,atom2) + ul2(3) = u_l(3,atom2) +#endif + + call getElementProperty(atypes, me1, "dipole_moment", mu1) + call getElementProperty(atypes, me2, "dipole_moment", mu2) + + +#ifdef IS_MPI + rf_Row(1,atom1) = rf_Row(1,atom1) + ul2(1)*mu2*taper + rf_Row(2,atom1) = rf_Row(2,atom1) + ul2(2)*mu2*taper + rf_Row(3,atom1) = rf_Row(3,atom1) + ul2(3)*mu2*taper + + rf_Col(1,atom2) = rf_Col(1,atom2) + ul1(1)*mu2*taper + rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu2*taper + rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu2*taper +#else + rf(1,atom1) = rf(1,atom1) + ul2(1)*mu2*taper + rf(2,atom1) = rf(2,atom1) + ul2(2)*mu2*taper + rf(3,atom1) = rf(3,atom1) + ul2(3)*mu2*taper + + rf(1,atom2) = rf(1,atom2) + ul1(1)*mu2*taper + rf(2,atom2) = rf(2,atom2) + ul1(2)*mu2*taper + rf(3,atom2) = rf(3,atom2) + ul1(3)*mu2*taper +#endif + + endif + return + end subroutine accumulate_rf -end subroutine reaction_field - -subroutine rf_correct_forces(atom1, atom2, dx, dy, dz, rij) - include 'sizes.inc' - include 'simulation.inc' - - integer atom1, atom2 - double precision dtdr, rrfsq, prerf, rij - double precision dudx, dudy, dudz, u1dotu2, dx, dy, dz - - rrfsq = rrf * rrf - prerf = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) + subroutine accumulate_self_rf(atom1, mu1, u_l) + + integer, intent(in) :: atom1 + real(kind=dp), intent(in) :: mu1 + real(kind=dp), dimension(:,:) :: u_l + + !! should work for both MPI and non-MPI version since this is not pairwise. + rf(1,atom1) = rf(1,atom1) + u_l(1,atom1)*mu1 + rf(2,atom1) = rf(2,atom1) + u_l(2,atom1)*mu1 + rf(3,atom1) = rf(3,atom1) + u_l(3,atom1)*mu1 + + return + end subroutine accumulate_self_rf - if (rij.le.rrf) then + subroutine reaction_field_final(a1, mu1, u_l, rfpot, t, do_pot) + + integer, intent(in) :: a1 + real (kind=dp), intent(in) :: mu1 + real (kind=dp), intent(inout) :: rfpot + logical, intent(in) :: do_pot + real (kind = dp), dimension(:,:) :: u_l + real (kind = dp), dimension(:,:) :: t - ! cubic taper - if (rij.lt.rt) then - dtdr = 0.0d0 - else - dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3) - endif + if (.not.rf_initialized) then + write(default_error,*) 'Reaction field not initialized!' + return + endif - u1dotu2 = ulx(atom1)*ulx(atom2) + uly(atom1)*uly(atom2) + & - ulz(atom1)*ulz(atom2) - - dudx = - prerf*mu(atom1)*mu(atom2)*u1dotu2*dtdr*dx/rij - dudy = - prerf*mu(atom1)*mu(atom2)*u1dotu2*dtdr*dy/rij - dudz = - prerf*mu(atom1)*mu(atom2)*u1dotu2*dtdr*dz/rij - - flx(atom1) = flx(atom1) + dudx - fly(atom1) = fly(atom1) + dudy - flz(atom1) = flz(atom1) + dudz + ! compute torques on dipoles: + ! pre converts from mu in units of debye to kcal/mol - flx(atom2) = flx(atom2) - dudx - fly(atom2) = fly(atom2) - dudy - flz(atom2) = flz(atom2) - dudz + ! The torque contribution is dipole cross reaction_field + + t(1,a1) = t(1,a1) + pre*mu1*(u_l(2,a1)*rf(3,a1) - u_l(3,a1)*rf(2,a1)) + t(2,a1) = t(2,a1) + pre*mu1*(u_l(3,a1)*rf(1,a1) - u_l(1,a1)*rf(3,a1)) + t(3,a1) = t(3,a1) + pre*mu1*(u_l(1,a1)*rf(2,a1) - u_l(2,a1)*rf(1,a1)) + + ! the potential contribution is -1/2 dipole dot reaction_field + + if (do_pot) then + rfpot = rfpot - 0.5d0 * pre * mu1 * & + (rf(1,a1)*u_l(1,a1) + rf(2,a1)*u_l(2,a1) + rf(3,a1)*u_l(3,a1)) + endif + + return + end subroutine reaction_field_final + + subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, f, do_stress) + + integer, intent(in) :: atom1, atom2 + real(kind=dp), dimension(3), intent(in) :: d + real(kind=dp), intent(in) :: rij + real( kind = dp ), dimension(:,:) :: u_l + real( kind = dp ), dimension(:,:) :: f + logical, intent(in) :: do_stress + + real (kind = dp), dimension(3) :: ul1 + real (kind = dp), dimension(3) :: ul2 + real (kind = dp) :: dtdr + real (kind = dp) :: dudx, dudy, dudz, u1dotu2 + integer :: me1, me2 + real (kind = dp) :: mu1, mu2 + + if (.not.rf_initialized) then + write(default_error,*) 'Reaction field not initialized!' + return + endif - ! add contribution to the virial as well - virial = virial + ( dx*dudx + dy*dudy + dz*dudz ) - - endif - - return -end subroutine rf_correct_forces + if (rij.le.rrf) then + + if (rij.lt.rt) then + dtdr = 0.0d0 + else + dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3) + endif + +#ifdef IS_MPI + me1 = atid_Row(atom1) + ul1(1) = u_l_Row(1,atom1) + ul1(2) = u_l_Row(2,atom1) + ul1(3) = u_l_Row(3,atom1) + + me2 = atid_Col(atom2) + ul2(1) = u_l_Col(1,atom2) + ul2(2) = u_l_Col(2,atom2) + ul2(3) = u_l_Col(3,atom2) +#else + me1 = atid(atom1) + ul1(1) = u_l(1,atom1) + ul1(2) = u_l(2,atom1) + ul1(3) = u_l(3,atom1) + + me2 = atid(atom2) + ul2(1) = u_l(1,atom2) + ul2(2) = u_l(2,atom2) + ul2(3) = u_l(3,atom2) +#endif + + call getElementProperty(atypes, me1, "dipole_moment", mu1) + call getElementProperty(atypes, me2, "dipole_moment", mu2) + + u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3) + + dudx = - pre*mu1*mu2*u1dotu2*dtdr*d(1)/rij + dudy = - pre*mu1*mu2*u1dotu2*dtdr*d(2)/rij + dudz = - pre*mu1*mu2*u1dotu2*dtdr*d(3)/rij + +#ifdef IS_MPI + f_Row(1,atom1) = f_Row(1,atom1) + dudx + f_Row(2,atom1) = f_Row(2,atom1) + dudy + f_Row(3,atom1) = f_Row(3,atom1) + dudz + + f_Col(1,atom2) = f_Col(1,atom2) - dudx + f_Col(2,atom2) = f_Col(2,atom2) - dudy + f_Col(3,atom2) = f_Col(3,atom2) - dudz +#else + f(1,atom1) = f(1,atom1) + dudx + f(2,atom1) = f(2,atom1) + dudy + f(3,atom1) = f(3,atom1) + dudz + + f(1,atom2) = f(1,atom2) - dudx + f(2,atom2) = f(2,atom2) - dudy + f(3,atom2) = f(3,atom2) - dudz +#endif + + if (do_stress) then + tau_Temp(1) = tau_Temp(1) + dudx * d(1) + tau_Temp(2) = tau_Temp(2) + dudx * d(2) + tau_Temp(3) = tau_Temp(3) + dudx * d(3) + tau_Temp(4) = tau_Temp(4) + dudy * d(1) + tau_Temp(5) = tau_Temp(5) + dudy * d(2) + tau_Temp(6) = tau_Temp(6) + dudy * d(3) + tau_Temp(7) = tau_Temp(7) + dudz * d(1) + tau_Temp(8) = tau_Temp(8) + dudz * d(2) + tau_Temp(9) = tau_Temp(9) + dudz * d(3) + virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) + endif + endif + + return + end subroutine rf_correct_forces +end module reaction_field