--- trunk/OOPSE/libmdtools/calc_LJ_FF.F90 2003/04/09 04:06:43 483 +++ trunk/OOPSE/libmdtools/calc_LJ_FF.F90 2004/06/01 21:45:22 1217 @@ -2,11 +2,12 @@ !! Corresponds to the force field defined in lj_FF.cpp !! @author Charles F. Vardeman II !! @author Matthew Meineke -!! @version $Id: calc_LJ_FF.F90,v 1.4 2003-04-09 04:06:43 gezelter Exp $, $Date: 2003-04-09 04:06:43 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $ +!! @version $Id: calc_LJ_FF.F90,v 1.22 2004-06-01 21:45:22 gezelter Exp $, $Date: 2004-06-01 21:45:22 $, $Name: not supported by cvs2svn $, $Revision: 1.22 $ module lj use definitions use atype_module + use switcheroo use vector_class use simulation #ifdef IS_MPI @@ -21,7 +22,10 @@ module lj #include "fForceField.h" integer, save :: LJ_Mixing_Policy - integer, save :: LJ_rcut + real(kind=DP), save :: LJ_rcut + logical, save :: havePolicy = .false. + logical, save :: haveCut = .false. + logical, save :: LJ_do_shift = .false. !! Logical has lj force field module been initialized? @@ -29,7 +33,7 @@ module lj !! Public methods and data public :: init_LJ_FF - public :: LJ_new_rcut + public :: setCutoffLJ public :: do_lj_pair type :: lj_mixed_params @@ -42,20 +46,17 @@ module lj !! real ( kind = dp ) :: tp6 real ( kind = dp ) :: tp12 - real ( kind = dp ) :: delta = 0.0_dp - end type lj_mixed_params type (lj_mixed_params), dimension(:,:), pointer :: ljMixed contains - subroutine init_LJ_FF(mix_Policy, rcut, status) + subroutine init_LJ_FF(mix_Policy, status) integer, intent(in) :: mix_Policy integer, intent(out) :: status integer :: myStatus - real(kind=dp) :: rcut if (mix_Policy == LB_MIXING_RULE) then LJ_Mixing_Policy = LB_MIXING_RULE @@ -68,34 +69,50 @@ contains return endif endif - - LJ_rcut = rcut - status = 0 - call createMixingList(myStatus) - if (myStatus /= 0) then - status = -1 - return + havePolicy = .true. + + if (haveCut) then + status = 0 + call createMixingList(myStatus) + if (myStatus /= 0) then + status = -1 + return + end if + + LJ_FF_initialized = .true. end if - - LJ_FF_initialized = .true. - + end subroutine init_LJ_FF - subroutine LJ_new_rcut(rcut, status) + subroutine setCutoffLJ(rcut, do_shift, status) + logical, intent(in):: do_shift integer :: status, myStatus real(kind=dp) :: rcut - LJ_rcut = rcut +#define __FORTRAN90 +#include "fSwitchingFunction.h" + status = 0 - call createMixingList(myStatus) - if (myStatus /= 0) then - status = -1 - return - end if + LJ_rcut = rcut + LJ_do_shift = do_shift + call set_switch(LJ_SWITCH, rcut, rcut) + haveCut = .true. + + if (havePolicy) then + status = 0 + call createMixingList(myStatus) + if (myStatus /= 0) then + status = -1 + return + end if + + LJ_FF_initialized = .true. + end if + return - end subroutine LJ_new_rcut + end subroutine setCutoffLJ subroutine createMixingList(status) integer :: nAtypes @@ -105,6 +122,7 @@ contains real ( kind = dp ) :: mySigma_i,mySigma_j real ( kind = dp ) :: myEpsilon_i,myEpsilon_j real ( kind = dp ) :: rcut6 + logical :: I_isLJ, J_isLJ status = 0 nAtypes = getSize(atypes) @@ -115,75 +133,87 @@ contains if (.not. associated(ljMixed)) then allocate(ljMixed(nAtypes, nAtypes)) - else - status = -1 - return - end if + endif rcut6 = LJ_rcut**6 +! This loops through all atypes, even those that don't support LJ forces. do i = 1, nAtypes - call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i) - call getElementProperty(atypes, i, "lj_sigma", mySigma_i) - ! do self mixing rule - ljMixed(i,i)%sigma = mySigma_i - - ljMixed(i,i)%sigma6 = (ljMixed(i,i)%sigma) ** 6 - ljMixed(i,i)%tp6 = ljMixed(i,i)%sigma6/rcut6 + call getElementProperty(atypes, i, "is_LJ", I_isLJ) - ljMixed(i,i)%tp12 = (ljMixed(i,i)%tp6) ** 2 + if (I_isLJ) then - ljMixed(i,i)%epsilon = myEpsilon_i - - ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * & - (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6) - - do j = i + 1, nAtypes - call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j) - call getElementProperty(atypes,j,"lj_sigma", mySigma_j) + call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i) + call getElementProperty(atypes, i, "lj_sigma", mySigma_i) + ! do self mixing rule + ljMixed(i,i)%sigma = mySigma_i - ljMixed(i,j)%sigma = & - calcLJMix("sigma",mySigma_i, & - mySigma_j) + ljMixed(i,i)%sigma6 = (ljMixed(i,i)%sigma) ** 6 - ljMixed(i,j)%sigma6 = & - (ljMixed(i,j)%sigma)**6 + ljMixed(i,i)%tp6 = (ljMixed(i,i)%sigma6)/rcut6 + ljMixed(i,i)%tp12 = (ljMixed(i,i)%tp6) ** 2 - ljMixed(i,j)%tp6 = ljMixed(i,j)%sigma6/rcut6 - ljMixed(i,j)%tp12 = (ljMixed(i,j)%tp6) ** 2 + ljMixed(i,i)%epsilon = myEpsilon_i + ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * & + (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6) - ljMixed(i,j)%epsilon = & - calcLJMix("epsilon",myEpsilon_i, & - myEpsilon_j) - - ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * & - (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6) - - - ljMixed(j,i)%sigma = ljMixed(i,j)%sigma - ljMixed(j,i)%sigma6 = ljMixed(i,j)%sigma6 - ljMixed(j,i)%tp6 = ljMixed(i,j)%tp6 - ljMixed(j,i)%tp12 = ljMixed(i,j)%tp12 - ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon - ljMixed(j,i)%delta = ljMixed(i,j)%delta - - end do + do j = i + 1, nAtypes + + call getElementProperty(atypes, j, "is_LJ", J_isLJ) + + if (J_isLJ) then + + call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j) + call getElementProperty(atypes,j,"lj_sigma", mySigma_j) + + ljMixed(i,j)%sigma = & + calcLJMix("sigma",mySigma_i, & + mySigma_j) + + ljMixed(i,j)%sigma6 = & + (ljMixed(i,j)%sigma)**6 + + + ljMixed(i,j)%tp6 = ljMixed(i,j)%sigma6/rcut6 + + ljMixed(i,j)%tp12 = (ljMixed(i,j)%tp6) ** 2 + + + ljMixed(i,j)%epsilon = & + calcLJMix("epsilon",myEpsilon_i, & + myEpsilon_j) + + ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * & + (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6) + + + ljMixed(j,i)%sigma = ljMixed(i,j)%sigma + ljMixed(j,i)%sigma6 = ljMixed(i,j)%sigma6 + ljMixed(j,i)%tp6 = ljMixed(i,j)%tp6 + ljMixed(j,i)%tp12 = ljMixed(i,j)%tp12 + ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon + ljMixed(j,i)%delta = ljMixed(i,j)%delta + endif + end do + endif end do end subroutine createMixingList - subroutine do_lj_pair(atom1, atom2, d, rij, r2, pot, f, do_pot, do_stress) + subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, & + pot, f, do_pot) integer, intent(in) :: atom1, atom2 real( kind = dp ), intent(in) :: rij, r2 - real( kind = dp ) :: pot - real( kind = dp ), dimension(3,getNlocal()) :: f + real( kind = dp ) :: pot, sw, vpair + real( kind = dp ), dimension(3,nLocal) :: f real( kind = dp ), intent(in), dimension(3) :: d - logical, intent(in) :: do_pot, do_stress + real( kind = dp ), intent(inout), dimension(3) :: fpair + logical, intent(in) :: do_pot ! local Variables real( kind = dp ) :: drdx, drdy, drdz @@ -195,90 +225,91 @@ contains real( kind = dp ) :: t6 real( kind = dp ) :: t12 real( kind = dp ) :: delta + integer :: id1, id2 - - if (rij.lt.LJ_rcut) then - - ! Look up the correct parameters in the mixing matrix + ! Look up the correct parameters in the mixing matrix #ifdef IS_MPI - sigma6 = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6 - epsilon = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon - delta = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta + sigma6 = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6 + epsilon = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon + delta = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta #else - sigma6 = ljMixed(atid(atom1),atid(atom2))%sigma6 - epsilon = ljMixed(atid(atom1),atid(atom2))%epsilon - delta = ljMixed(atid(atom1),atid(atom2))%delta + sigma6 = ljMixed(atid(atom1),atid(atom2))%sigma6 + epsilon = ljMixed(atid(atom1),atid(atom2))%epsilon + delta = ljMixed(atid(atom1),atid(atom2))%delta #endif + + r6 = r2 * r2 * r2 + + t6 = sigma6/ r6 + t12 = t6 * t6 + + pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + if (LJ_do_shift) then + pot_temp = pot_temp + delta + endif + + vpair = vpair + pot_temp - r6 = r2 * r2 * r2 + dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij - t6 = sigma6/ r6 - t12 = t6 * t6 - - pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta + drdx = d(1) / rij + drdy = d(2) / rij + drdz = d(3) / rij - dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij + fx = dudr * drdx + fy = dudr * drdy + fz = dudr * drdz + - drdx = d(1) / rij - drdy = d(2) / rij - drdz = d(3) / rij - - fx = dudr * drdx - fy = dudr * drdy - fz = dudr * drdz - #ifdef IS_MPI - if (do_pot) then - pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5 - pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5 - endif - - f_Row(1,atom1) = f_Row(1,atom1) + fx - f_Row(2,atom1) = f_Row(2,atom1) + fy - f_Row(3,atom1) = f_Row(3,atom1) + fz - - f_Col(1,atom2) = f_Col(1,atom2) - fx - f_Col(2,atom2) = f_Col(2,atom2) - fy - f_Col(3,atom2) = f_Col(3,atom2) - fz - + if (do_pot) then + pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5 + pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5 + endif + + f_Row(1,atom1) = f_Row(1,atom1) + fx + f_Row(2,atom1) = f_Row(2,atom1) + fy + f_Row(3,atom1) = f_Row(3,atom1) + fz + + f_Col(1,atom2) = f_Col(1,atom2) - fx + f_Col(2,atom2) = f_Col(2,atom2) - fy + f_Col(3,atom2) = f_Col(3,atom2) - fz + #else - if (do_pot) pot = pot + pot_temp + if (do_pot) pot = pot + sw*pot_temp - f(1,atom1) = f(1,atom1) + fx - f(2,atom1) = f(2,atom1) + fy - f(3,atom1) = f(3,atom1) + fz - - f(1,atom2) = f(1,atom2) - fx - f(2,atom2) = f(2,atom2) - fy - f(3,atom2) = f(3,atom2) - fz + f(1,atom1) = f(1,atom1) + fx + f(2,atom1) = f(2,atom1) + fy + f(3,atom1) = f(3,atom1) + fz + + f(1,atom2) = f(1,atom2) - fx + f(2,atom2) = f(2,atom2) - fy + f(3,atom2) = f(3,atom2) - fz #endif - - if (do_stress) then + +#ifdef IS_MPI + id1 = AtomRowToGlobal(atom1) + id2 = AtomColToGlobal(atom2) +#else + id1 = atom1 + id2 = atom2 +#endif - if (molMembershipList(atom1) .ne. molMembershipList(atom2)) then - tau_Temp(1) = tau_Temp(1) + fx * d(1) - tau_Temp(2) = tau_Temp(2) + fx * d(2) - tau_Temp(3) = tau_Temp(3) + fx * d(3) - tau_Temp(4) = tau_Temp(4) + fy * d(1) - tau_Temp(5) = tau_Temp(5) + fy * d(2) - tau_Temp(6) = tau_Temp(6) + fy * d(3) - tau_Temp(7) = tau_Temp(7) + fz * d(1) - tau_Temp(8) = tau_Temp(8) + fz * d(2) - tau_Temp(9) = tau_Temp(9) + fz * d(3) - virial_Temp = virial_Temp + & - (tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) - endif - - endif + if (molMembershipList(id1) .ne. molMembershipList(id2)) then + fpair(1) = fpair(1) + fx + fpair(2) = fpair(2) + fy + fpair(3) = fpair(3) + fz + endif + return - + end subroutine do_lj_pair - - + + !! Calculates the mixing for sigma or epslon - + function calcLJMix(thisParam,param1,param2,status) result(myMixParam) character(len=*) :: thisParam real(kind = dp) :: param1 @@ -291,7 +322,7 @@ contains if (present(status)) status = 0 select case (LJ_Mixing_Policy) - case (LB_MIXING_RULE) + case (1) select case (thisParam) case ("sigma") myMixParam = 0.5_dp * (param1 + param2)