--- trunk/OOPSE/libmdtools/calc_LJ_FF.F90 2004/01/05 22:49:14 898 +++ trunk/OOPSE/libmdtools/calc_LJ_FF.F90 2004/05/27 00:48:12 1198 @@ -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.14 2004-01-05 22:49:14 chuckv Exp $, $Date: 2004-01-05 22:49:14 $, $Name: not supported by cvs2svn $, $Revision: 1.14 $ +!! @version $Id: calc_LJ_FF.F90,v 1.21 2004-05-27 00:48:12 tim Exp $, $Date: 2004-05-27 00:48:12 $, $Name: not supported by cvs2svn $, $Revision: 1.21 $ module lj use definitions use atype_module + use switcheroo use vector_class use simulation #ifdef IS_MPI @@ -22,8 +23,9 @@ module lj integer, save :: LJ_Mixing_Policy real(kind=DP), save :: LJ_rcut - logical, save :: havePolicy = .false., haveCut = .false. - + logical, save :: havePolicy = .false. + logical, save :: haveCut = .false. + logical, save :: LJ_do_shift = .false. !! Logical has lj force field module been initialized? @@ -44,10 +46,7 @@ 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 @@ -86,13 +85,19 @@ contains end subroutine init_LJ_FF - subroutine setCutoffLJ(rcut, status) + subroutine setCutoffLJ(rcut, do_shift, status) + logical, intent(in):: do_shift integer :: status, myStatus real(kind=dp) :: rcut - + +#define __FORTRAN90 +#include "fSwitchingFunction.h" + status = 0 LJ_rcut = rcut + LJ_do_shift = do_shift + call set_switch(LJ_SWITCH, rcut, rcut) haveCut = .true. if (havePolicy) then @@ -131,6 +136,7 @@ contains 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) @@ -187,14 +193,16 @@ contains 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 ) :: 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 @@ -208,103 +216,89 @@ contains 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 = 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 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 + 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