--- trunk/OOPSE/libmdtools/calc_reaction_field.F90 2004/01/05 22:49:14 898 +++ trunk/OOPSE/libmdtools/calc_reaction_field.F90 2004/05/24 21:03:30 1192 @@ -4,6 +4,7 @@ module reaction_field use atype_module use vector_class use simulation + use status #ifdef IS_MPI use mpiSimulation #endif @@ -16,9 +17,16 @@ module reaction_field real(kind=dp), save :: dielect = 1.0_dp real(kind=dp), save :: rrfsq = 1.0_dp real(kind=dp), save :: pre - logical, save :: rf_initialized = .false., haveCuts = .false. - logical, save :: haveDie = .false. + logical, save :: haveCutoffs = .false. + logical, save :: haveMomentMap = .false. + logical, save :: haveDielectric = .false. + type :: MomentList + real(kind=DP) :: dipole_moment = 0.0_DP + end type MomentList + + type(MomentList), dimension(:),allocatable :: MomentMap + PUBLIC::initialize_rf PUBLIC::setCutoffsRF PUBLIC::accumulate_rf @@ -35,8 +43,7 @@ contains pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) - haveDie = .true. - if (haveCuts) rf_initialized = .true. + haveDielectric = .true. return end subroutine initialize_rf @@ -51,82 +58,118 @@ contains rrfsq = rrf * rrf pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) - haveCuts = .true. - if (haveDie) rf_initialized = .true. + haveCutoffs = .true. end subroutine setCutoffsRF - - subroutine accumulate_rf(atom1, atom2, rij, u_l) + subroutine createMomentMap(status) + integer :: nAtypes + integer :: status + integer :: i + real (kind=DP) :: thisDP + logical :: thisProperty + status = 0 + + nAtypes = getSize(atypes) + + if (nAtypes == 0) then + status = -1 + return + end if + + if (.not. allocated(MomentMap)) then + allocate(MomentMap(nAtypes)) + endif + + do i = 1, nAtypes + + call getElementProperty(atypes, i, "is_DP", thisProperty) + + if (thisProperty) then + call getElementProperty(atypes, i, "dipole_moment", thisDP) + MomentMap(i)%dipole_moment = thisDP + endif + + end do + + haveMomentMap = .true. + + end subroutine createMomentMap + + 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 - if (.not.rf_initialized) then + integer :: localError + + if ((.not.haveDielectric).or.(.not.haveCutoffs)) then write(default_error,*) 'Reaction field not initialized!' return endif + + if (.not.haveMomentMap) then + localError = 0 + call createMomentMap(localError) + if ( localError .ne. 0 ) then + call handleError("reaction-field", "MomentMap creation failed!") + return + 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 - - call getElementProperty(atypes, me1, "dipole_moment", mu1) - call getElementProperty(atypes, me2, "dipole_moment", mu2) - - + + 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 + 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(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 @@ -153,11 +196,22 @@ contains real (kind = dp), dimension(3,nLocal) :: u_l real (kind = dp), dimension(3,nLocal) :: t - if (.not.rf_initialized) then + integer :: localError + + if ((.not.haveDielectric).or.(.not.haveCutoffs)) then write(default_error,*) 'Reaction field not initialized!' return endif + if (.not.haveMomentMap) then + localError = 0 + call createMomentMap(localError) + if ( localError .ne. 0 ) then + call handleError("reaction-field", "MomentMap creation failed!") + return + end if + endif + ! compute torques on dipoles: ! pre converts from mu in units of debye to kcal/mol @@ -173,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 @@ -193,17 +247,28 @@ contains integer :: me1, me2, id1, id2 real (kind = dp) :: mu1, mu2 - if (.not.rf_initialized) then + integer :: localError + + if ((.not.haveDielectric).or.(.not.haveCutoffs)) then write(default_error,*) 'Reaction field not initialized!' return endif + if (.not.haveMomentMap) then + localError = 0 + call createMomentMap(localError) + if ( localError .ne. 0 ) then + call handleError("reaction-field", "MomentMap creation failed!") + return + end if + endif + if (rij.le.rrf) then 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 @@ -229,8 +294,8 @@ contains ul2(3) = u_l(3,atom2) #endif - call getElementProperty(atypes, me1, "dipole_moment", mu1) - call getElementProperty(atypes, me2, "dipole_moment", mu2) + mu1 = MomentMap(me1)%dipole_moment + mu2 = MomentMap(me2)%dipole_moment u1dotu2 = ul1(1)*ul2(1) + ul1(2)*ul2(2) + ul1(3)*ul2(3) @@ -255,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 = tagRow(atom1) + id2 = tagColumn(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