--- branches/mmeineke/OOPSE/libmdtools/calc_LJ_FF.F90 2003/03/21 17:42:12 377 +++ trunk/OOPSE/libmdtools/calc_LJ_FF.F90 2004/02/09 14:48:57 1041 @@ -2,12 +2,13 @@ !! 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.1.1.1 2003-03-21 17:42:12 mmeineke Exp $, $Date: 2003-03-21 17:42:12 $, $Name: not supported by cvs2svn $, $Revision: 1.1.1.1 $ +!! @version $Id: calc_LJ_FF.F90,v 1.17 2004-02-09 14:48:57 chrisfen Exp $, $Date: 2004-02-09 14:48:57 $, $Name: not supported by cvs2svn $, $Revision: 1.17 $ module lj use definitions use atype_module use vector_class + use simulation #ifdef IS_MPI use mpiSimulation #endif @@ -20,7 +21,9 @@ module lj #include "fForceField.h" integer, save :: LJ_Mixing_Policy - integer, save :: LJ_rcut + real(kind=DP), save :: LJ_rcut + logical, save :: havePolicy = .false., haveCut = .false. + !! Logical has lj force field module been initialized? @@ -28,7 +31,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 @@ -44,17 +47,17 @@ module lj 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 @@ -67,34 +70,46 @@ 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, status) integer :: status, myStatus real(kind=dp) :: rcut + + status = 0 LJ_rcut = rcut - status = 0 - call createMixingList(myStatus) - if (myStatus /= 0) then - status = -1 - return - end if +!!$ ! ATTENTION! This is a hardwiring of rcut! +!!$ LJ_rcut = 9.0d0 + 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 @@ -114,13 +129,11 @@ 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) @@ -128,11 +141,13 @@ contains ! 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 + ljMixed(i,i)%sigma6 = (ljMixed(i,i)%sigma) ** 6 + ljMixed(i,i)%tp6 = (ljMixed(i,i)%sigma6)/rcut6 + ljMixed(i,i)%tp12 = (ljMixed(i,i)%tp6) ** 2 + ljMixed(i,i)%epsilon = myEpsilon_i ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * & @@ -180,7 +195,7 @@ contains integer, intent(in) :: atom1, atom2 real( kind = dp ), intent(in) :: rij, r2 real( kind = dp ) :: pot - real( kind = dp ), dimension(:,:) :: f + real( kind = dp ), dimension(3,nLocal) :: f real( kind = dp ), intent(in), dimension(3) :: d logical, intent(in) :: do_pot, do_stress @@ -194,8 +209,8 @@ 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 @@ -218,9 +233,9 @@ contains dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij - drdx = -d(1) / rij - drdy = -d(2) / rij - drdz = -d(3) / rij + drdx = d(1) / rij + drdy = d(2) / rij + drdz = d(3) / rij fx = dudr * drdx fy = dudr * drdy @@ -253,16 +268,35 @@ contains #endif if (do_stress) 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)) + +#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 fx, fy, fz are the force on atom i, we need a + ! negative sign here: + + tau_Temp(1) = tau_Temp(1) - d(1) * fx + tau_Temp(2) = tau_Temp(2) - d(1) * fy + tau_Temp(3) = tau_Temp(3) - d(1) * fz + tau_Temp(4) = tau_Temp(4) - d(2) * fx + tau_Temp(5) = tau_Temp(5) - d(2) * fy + tau_Temp(6) = tau_Temp(6) - d(2) * fz + tau_Temp(7) = tau_Temp(7) - d(3) * fx + tau_Temp(8) = tau_Temp(8) - d(3) * fy + tau_Temp(9) = tau_Temp(9) - d(3) * fz + + virial_Temp = virial_Temp + & + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) + + endif endif endif @@ -285,7 +319,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)