--- trunk/OOPSE/libmdtools/calc_reaction_field.F90 2003/04/04 22:22:30 460 +++ trunk/OOPSE/libmdtools/calc_reaction_field.F90 2004/01/30 15:01:09 999 @@ -4,61 +4,132 @@ module reaction_field use atype_module use vector_class use simulation + use status #ifdef IS_MPI use mpiSimulation #endif implicit none - real(kind=dp), save :: rrf + PRIVATE + + real(kind=dp), save :: rrf = 1.0_dp real(kind=dp), save :: rt - real(kind=dp), save :: dielect - real(kind=dp), save :: rrfsq + 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. + 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 + PUBLIC::accumulate_self_rf + PUBLIC::reaction_field_final + PUBLIC::rf_correct_forces + contains - subroutine initialize_rf(this_rrf, this_rt, this_dielect) - real(kind=dp), intent(in) :: this_rrf, this_rt, this_dielect - - rrf = this_rrf - rt = this_rt + subroutine initialize_rf(this_dielect) + real(kind=dp), intent(in) :: this_dielect + dielect = this_dielect + + pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) + haveDielectric = .true. + + return + end subroutine initialize_rf + + subroutine setCutoffsRF( this_rrf, this_rt ) + + real(kind=dp), intent(in) :: this_rrf, this_rt + + rrf = this_rrf + rt = this_rt + rrfsq = rrf * rrf pre = 14.38362d0*2.0d0*(dielect-1.0d0)/((2.0d0*dielect+1.0d0)*rrfsq*rrf) + haveCutoffs = .true. - write(*,*) 'rrf = ', rrf - write(*,*) 'rt = ', rt - write(*,*) 'dielect = ', dielect - write(*,*) 'pre = ', pre - rf_initialized = .true. + end subroutine setCutoffsRF + + subroutine createMomentMap(status) + integer :: nAtypes + integer :: status + integer :: i + real (kind=DP) :: thisDP + logical :: thisProperty + + status = 0 + + nAtypes = getSize(atypes) - return - end subroutine initialize_rf - + 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) integer, intent(in) :: atom1, atom2 real (kind = dp), intent(in) :: rij - real (kind = dp), dimension(3,getNlocal()) :: u_l + real (kind = dp), dimension(3,nLocal) :: u_l integer :: me1, me2 real (kind = dp) :: taper, 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 @@ -84,10 +155,9 @@ 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 - #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 @@ -114,7 +184,7 @@ contains integer, intent(in) :: atom1 real(kind=dp), intent(in) :: mu1 - real(kind=dp), dimension(3,getNlocal()) :: u_l + real(kind=dp), dimension(3,nLocal) :: 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 @@ -130,14 +200,25 @@ contains real (kind=dp), intent(in) :: mu1 real (kind=dp), intent(inout) :: rfpot logical, intent(in) :: do_pot - real (kind = dp), dimension(3,getNlocal()) :: u_l - real (kind = dp), dimension(3,getNlocal()) :: t + 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 @@ -153,7 +234,7 @@ 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 @@ -162,27 +243,39 @@ contains integer, intent(in) :: atom1, atom2 real(kind=dp), dimension(3), intent(in) :: d real(kind=dp), intent(in) :: rij - real( kind = dp ), dimension(3,getNlocal()) :: u_l - real( kind = dp ), dimension(3,getNlocal()) :: f + 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) :: ul1 real (kind = dp), dimension(3) :: ul2 real (kind = dp) :: dtdr real (kind = dp) :: dudx, dudy, dudz, u1dotu2 - integer :: me1, me2 + 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' dtdr = 6.0d0*(rij*rij - rij*rt - rij*rrf +rrf*rt)/((rrf-rt)**3) endif @@ -208,8 +301,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) @@ -236,16 +329,34 @@ contains #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)) + +#ifdef IS_MPI + id1 = tagRow(atom1) + id2 = tagColumn(atom2) +#else + 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