--- trunk/OOPSE/libmdtools/calc_reaction_field.F90 2004/01/06 19:49:18 901 +++ trunk/OOPSE/libmdtools/calc_reaction_field.F90 2004/05/27 00:48:12 1198 @@ -97,14 +97,15 @@ contains end subroutine createMomentMap - subroutine accumulate_rf(atom1, atom2, rij, u_l) + subroutine accumulate_rf(atom1, atom2, rij, u_l, taper) integer, intent(in) :: atom1, atom2 real (kind = dp), intent(in) :: rij real (kind = dp), dimension(3,nLocal) :: u_l integer :: me1, me2 - real (kind = dp) :: taper, mu1, mu2 + real (kind = dp), intent(in) :: taper + real (kind = dp):: mu1, mu2 real (kind = dp), dimension(3) :: ul1 real (kind = dp), dimension(3) :: ul2 @@ -124,59 +125,51 @@ contains end if endif - if (rij.le.rrf) then - - if (rij.lt.rt) then - taper = 1.0d0 - else - write(*,*) 'rf in taper region' - 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) + 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) + 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 - - mu1 = MomentMap(me1)%dipole_moment - mu2 = MomentMap(me2)%dipole_moment - + + mu1 = MomentMap(me1)%dipole_moment + mu2 = MomentMap(me2)%dipole_moment + #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)*mu1*taper - rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper - rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*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)*mu1*taper - rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper - rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper + 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)*mu1*taper + rf_Col(2,atom2) = rf_Col(2,atom2) + ul1(2)*mu1*taper + rf_Col(3,atom2) = rf_Col(3,atom2) + ul1(3)*mu1*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)*mu1*taper + rf(2,atom2) = rf(2,atom2) + ul1(2)*mu1*taper + rf(3,atom2) = rf(3,atom2) + ul1(3)*mu1*taper #endif - - endif + + return end subroutine accumulate_rf @@ -234,18 +227,18 @@ contains 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) + subroutine rf_correct_forces(atom1, atom2, d, rij, u_l, taper, f, fpair) integer, intent(in) :: atom1, atom2 real(kind=dp), dimension(3), intent(in) :: d - real(kind=dp), intent(in) :: rij + real(kind=dp), intent(in) :: rij, taper real( kind = dp ), dimension(3,nLocal) :: u_l real( kind = dp ), dimension(3,nLocal) :: f - logical, intent(in) :: do_stress + real( kind = dp ), dimension(3), intent(inout) :: fpair real (kind = dp), dimension(3) :: ul1 real (kind = dp), dimension(3) :: ul2 @@ -275,7 +268,7 @@ contains if (rij.lt.rt) then dtdr = 0.0d0 else - write(*,*) 'rf correct in taper region' + ! write(*,*) 'rf correct in taper region' dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3) endif @@ -327,39 +320,24 @@ contains f(2,atom2) = f(2,atom2) - dudy f(3,atom2) = f(3,atom2) - dudz #endif - - if (do_stress) then #ifdef IS_MPI - id1 = tagRow(atom1) - id2 = tagColumn(atom2) + id1 = AtomRowToGlobal(atom1) + id2 = AtomColToGlobal(atom2) #else - id1 = atom1 - id2 = atom2 + id1 = atom1 + id2 = atom2 #endif - - if (molMembershipList(id1) .ne. molMembershipList(id2)) then - - ! because the d vector is the rj - ri vector, and - ! because dudx, dudy, and dudz are the - ! (positive) force on atom i (negative on atom j) we need - ! a negative sign here: - - tau_Temp(1) = tau_Temp(1) - d(1) * dudx - tau_Temp(2) = tau_Temp(2) - d(1) * dudy - tau_Temp(3) = tau_Temp(3) - d(1) * dudz - tau_Temp(4) = tau_Temp(4) - d(2) * dudx - tau_Temp(5) = tau_Temp(5) - d(2) * dudy - tau_Temp(6) = tau_Temp(6) - d(2) * dudz - tau_Temp(7) = tau_Temp(7) - d(3) * dudx - tau_Temp(8) = tau_Temp(8) - d(3) * dudy - tau_Temp(9) = tau_Temp(9) - d(3) * dudz - virial_Temp = virial_Temp + & - (tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) - endif - endif - endif - + + if (molMembershipList(id1) .ne. molMembershipList(id2)) then + + fpair(1) = fpair(1) + dudx + fpair(2) = fpair(2) + dudy + fpair(3) = fpair(3) + dudz + + endif + + end if return end subroutine rf_correct_forces end module reaction_field