OOPSE/libmdtools/calc_LJ_FF.F90

!! Calculates Long Range forces Lennard-Jones interactions.
!! 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.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
  use mpiSimulation
#endif
  use force_globals

  implicit none
  PRIVATE

#define __FORTRAN90
#include "fForceField.h"

  integer, save :: LJ_Mixing_Policy
  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?
  
  logical, save :: LJ_FF_initialized = .false.
  
  !! Public methods and data
  public :: init_LJ_FF
  public :: setCutoffLJ
  public :: do_lj_pair
  
  type :: lj_mixed_params
     !! Lennard-Jones epsilon
     real ( kind = dp )  :: epsilon = 0.0_dp
     !! Lennard-Jones Sigma
     real ( kind = dp )  :: sigma = 0.0_dp
     !! Lennard-Jones Sigma to sixth
     real ( kind = dp )  :: sigma6 = 0.0_dp
     !! 
     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, status)
    integer, intent(in) :: mix_Policy
    integer, intent(out) :: status
    integer :: myStatus
    
    if (mix_Policy == LB_MIXING_RULE) then
       LJ_Mixing_Policy = LB_MIXING_RULE
    else
       if (mix_Policy == EXPLICIT_MIXING_RULE) then
          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
       else
          write(*,*) 'Unknown Mixing Policy!'
          status = -1
          return
       endif
    endif

    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
  
  end subroutine init_LJ_FF
  
  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
       status = 0
       call createMixingList(myStatus)
       if (myStatus /= 0) then
          status = -1
          return
       end if
       
       LJ_FF_initialized = .true.
    end if    
    
    return
  end subroutine setCutoffLJ
  
  subroutine createMixingList(status)
    integer :: nAtypes
    integer :: status
    integer :: i
    integer :: j
    real ( kind = dp ) :: mySigma_i,mySigma_j
    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
    real ( kind = dp ) :: rcut6
    status = 0
    
    nAtypes = getSize(atypes)
    if (nAtypes == 0) then
       status = -1
       return
    end if
        
    if (.not. associated(ljMixed)) then
       allocate(ljMixed(nAtypes, nAtypes))
    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

       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 * &
            (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)
          
          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
          
       end do
    end do
    
  end subroutine createMixingList
  
  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, sw, vpair
    real( kind = dp ), dimension(3,nLocal) :: f    
    real( kind = dp ), intent(in), dimension(3) :: d
    real( kind = dp ), intent(inout), dimension(3) :: fpair
    logical, intent(in) :: do_pot

    ! local Variables
    real( kind = dp ) :: drdx, drdy, drdz
    real( kind = dp ) :: fx, fy, fz
    real( kind = dp ) :: pot_temp, dudr
    real( kind = dp ) :: sigma6
    real( kind = dp ) :: epsilon
    real( kind = dp ) :: r6
    real( kind = dp ) :: t6
    real( kind = dp ) :: t12
    real( kind = dp ) :: delta
    integer :: id1, id2

    ! 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
#else
    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
       
    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
       
    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) + 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 + 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
#endif
        
#ifdef IS_MPI
    id1 = AtomRowToGlobal(atom1)
    id2 = AtomColToGlobal(atom2)
#else
    id1 = atom1
    id2 = atom2
#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
    real(kind = dp)  :: param2
    real(kind = dp ) :: myMixParam

    integer, optional :: status   

    myMixParam = 0.0_dp
    
    if (present(status)) status = 0
    select case (LJ_Mixing_Policy)
    case (1)
       select case (thisParam)
       case ("sigma")
          myMixParam = 0.5_dp * (param1 + param2)
       case ("epsilon")
          myMixParam = sqrt(param1 * param2)
       case default
          status = -1
       end select
    case default
       status = -1
    end select
  end function calcLJMix
  
end module lj
Revision:	1198
Committed:	Thu May 27 00:48:12 2004 UTC (20 years, 1 month ago) by tim
File size:	8551 byte(s)
Log Message:	in the progress of fixing MPI version of cutoff group
#	User	Rev	Content
1	mmeineke	377	!! Calculates Long Range forces Lennard-Jones interactions.
2			!! Corresponds to the force field defined in lj_FF.cpp
3			!! @author Charles F. Vardeman II
4			!! @author Matthew Meineke
5	tim	1198	!! @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 $
6	mmeineke	377
7			module lj
8			use definitions
9			use atype_module
10	gezelter	1150	use switcheroo
11	mmeineke	377	use vector_class
12	chuckv	460	use simulation
13	mmeineke	377	#ifdef IS_MPI
14			use mpiSimulation
15			#endif
16			use force_globals
17
18			implicit none
19			PRIVATE
20
21			#define __FORTRAN90
22			#include "fForceField.h"
23
24			integer, save :: LJ_Mixing_Policy
25	mmeineke	619	real(kind=DP), save :: LJ_rcut
26	gezelter	1160	logical, save :: havePolicy = .false.
27			logical, save :: haveCut = .false.
28			logical, save :: LJ_do_shift = .false.
29	mmeineke	377
30			!! Logical has lj force field module been initialized?
31
32			logical, save :: LJ_FF_initialized = .false.
33
34			!! Public methods and data
35			public :: init_LJ_FF
36	mmeineke	626	public :: setCutoffLJ
37	mmeineke	377	public :: do_lj_pair
38
39			type :: lj_mixed_params
40			!! Lennard-Jones epsilon
41			real ( kind = dp ) :: epsilon = 0.0_dp
42			!! Lennard-Jones Sigma
43			real ( kind = dp ) :: sigma = 0.0_dp
44			!! Lennard-Jones Sigma to sixth
45			real ( kind = dp ) :: sigma6 = 0.0_dp
46			!!
47			real ( kind = dp ) :: tp6
48			real ( kind = dp ) :: tp12
49			real ( kind = dp ) :: delta = 0.0_dp
50			end type lj_mixed_params
51
52			type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
53
54			contains
55
56	mmeineke	626	subroutine init_LJ_FF(mix_Policy, status)
57	mmeineke	377	integer, intent(in) :: mix_Policy
58			integer, intent(out) :: status
59			integer :: myStatus
60
61	gezelter	834	if (mix_Policy == LB_MIXING_RULE) then
62			LJ_Mixing_Policy = LB_MIXING_RULE
63	mmeineke	377	else
64	gezelter	834	if (mix_Policy == EXPLICIT_MIXING_RULE) then
65			LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
66	mmeineke	377	else
67			write(,) 'Unknown Mixing Policy!'
68			status = -1
69			return
70			endif
71			endif
72
73	mmeineke	626	havePolicy = .true.
74
75			if (haveCut) then
76			status = 0
77			call createMixingList(myStatus)
78			if (myStatus /= 0) then
79			status = -1
80			return
81			end if
82
83			LJ_FF_initialized = .true.
84	mmeineke	377	end if
85	mmeineke	626
86	mmeineke	377	end subroutine init_LJ_FF
87
88	gezelter	1160	subroutine setCutoffLJ(rcut, do_shift, status)
89			logical, intent(in):: do_shift
90	mmeineke	377	integer :: status, myStatus
91			real(kind=dp) :: rcut
92	gezelter	1150
93			#define __FORTRAN90
94			#include "fSwitchingFunction.h"
95
96	mmeineke	377	status = 0
97
98	chrisfen	1041	LJ_rcut = rcut
99	gezelter	1160	LJ_do_shift = do_shift
100	gezelter	1150	call set_switch(LJ_SWITCH, rcut, rcut)
101	mmeineke	626	haveCut = .true.
102
103			if (havePolicy) then
104			status = 0
105	mmeineke	619	call createMixingList(myStatus)
106			if (myStatus /= 0) then
107			status = -1
108			return
109			end if
110	mmeineke	626
111			LJ_FF_initialized = .true.
112			end if
113	mmeineke	619
114	mmeineke	377	return
115	mmeineke	626	end subroutine setCutoffLJ
116	mmeineke	377
117			subroutine createMixingList(status)
118			integer :: nAtypes
119			integer :: status
120			integer :: i
121			integer :: j
122			real ( kind = dp ) :: mySigma_i,mySigma_j
123			real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
124			real ( kind = dp ) :: rcut6
125			status = 0
126
127			nAtypes = getSize(atypes)
128			if (nAtypes == 0) then
129			status = -1
130			return
131			end if
132
133			if (.not. associated(ljMixed)) then
134			allocate(ljMixed(nAtypes, nAtypes))
135	chuckv	623	endif
136	mmeineke	377
137			rcut6 = LJ_rcut**6
138
139	chuckv	900	! This loops through all atypes, even those that don't support LJ forces.
140	mmeineke	377	do i = 1, nAtypes
141
142			call getElementProperty(atypes, i, "lj_epsilon", myEpsilon_i)
143			call getElementProperty(atypes, i, "lj_sigma", mySigma_i)
144			! do self mixing rule
145			ljMixed(i,i)%sigma = mySigma_i
146
147	mmeineke	790	ljMixed(i,i)%sigma6 = (ljMixed(i,i)%sigma) ** 6
148	mmeineke	377
149	mmeineke	790	ljMixed(i,i)%tp6 = (ljMixed(i,i)%sigma6)/rcut6
150
151	mmeineke	377	ljMixed(i,i)%tp12 = (ljMixed(i,i)%tp6) ** 2
152
153	mmeineke	790
154	mmeineke	377	ljMixed(i,i)%epsilon = myEpsilon_i
155
156			ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
157			(ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
158
159			do j = i + 1, nAtypes
160			call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
161			call getElementProperty(atypes,j,"lj_sigma", mySigma_j)
162
163			ljMixed(i,j)%sigma = &
164			calcLJMix("sigma",mySigma_i, &
165			mySigma_j)
166
167			ljMixed(i,j)%sigma6 = &
168			(ljMixed(i,j)%sigma)**6
169
170
171			ljMixed(i,j)%tp6 = ljMixed(i,j)%sigma6/rcut6
172
173			ljMixed(i,j)%tp12 = (ljMixed(i,j)%tp6) ** 2
174
175
176			ljMixed(i,j)%epsilon = &
177			calcLJMix("epsilon",myEpsilon_i, &
178			myEpsilon_j)
179
180			ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
181			(ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
182
183
184			ljMixed(j,i)%sigma = ljMixed(i,j)%sigma
185			ljMixed(j,i)%sigma6 = ljMixed(i,j)%sigma6
186			ljMixed(j,i)%tp6 = ljMixed(i,j)%tp6
187			ljMixed(j,i)%tp12 = ljMixed(i,j)%tp12
188			ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
189			ljMixed(j,i)%delta = ljMixed(i,j)%delta
190
191			end do
192			end do
193
194			end subroutine createMixingList
195
196	gezelter	1192	subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
197			pot, f, do_pot)
198	mmeineke	377
199			integer, intent(in) :: atom1, atom2
200			real( kind = dp ), intent(in) :: rij, r2
201	gezelter	1150	real( kind = dp ) :: pot, sw, vpair
202	chuckv	898	real( kind = dp ), dimension(3,nLocal) :: f
203	mmeineke	377	real( kind = dp ), intent(in), dimension(3) :: d
204	gezelter	1192	real( kind = dp ), intent(inout), dimension(3) :: fpair
205			logical, intent(in) :: do_pot
206	mmeineke	377
207			! local Variables
208			real( kind = dp ) :: drdx, drdy, drdz
209			real( kind = dp ) :: fx, fy, fz
210			real( kind = dp ) :: pot_temp, dudr
211			real( kind = dp ) :: sigma6
212			real( kind = dp ) :: epsilon
213			real( kind = dp ) :: r6
214			real( kind = dp ) :: t6
215			real( kind = dp ) :: t12
216			real( kind = dp ) :: delta
217	gezelter	490	integer :: id1, id2
218	mmeineke	377
219	gezelter	1160	! Look up the correct parameters in the mixing matrix
220	mmeineke	377	#ifdef IS_MPI
221	gezelter	1160	sigma6 = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
222			epsilon = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
223			delta = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
224	mmeineke	377	#else
225	gezelter	1160	sigma6 = ljMixed(atid(atom1),atid(atom2))%sigma6
226			epsilon = ljMixed(atid(atom1),atid(atom2))%epsilon
227			delta = ljMixed(atid(atom1),atid(atom2))%delta
228	mmeineke	377	#endif
229	gezelter	1160
230			r6 = r2 * r2 * r2
231
232			t6 = sigma6/ r6
233			t12 = t6 * t6
234
235			pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
236			if (LJ_do_shift) then
237			pot_temp = pot_temp + delta
238			endif
239
240			vpair = vpair + pot_temp
241	mmeineke	377
242	gezelter	1160	dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
243	mmeineke	377
244	gezelter	1160	drdx = d(1) / rij
245			drdy = d(2) / rij
246			drdz = d(3) / rij
247	mmeineke	377
248	gezelter	1160	fx = dudr * drdx
249			fy = dudr * drdy
250			fz = dudr * drdz
251
252	mmeineke	377
253			#ifdef IS_MPI
254	gezelter	1160	if (do_pot) then
255			pot_Row(atom1) = pot_Row(atom1) + swpot_temp0.5
256			pot_Col(atom2) = pot_Col(atom2) + swpot_temp0.5
257			endif
258
259			f_Row(1,atom1) = f_Row(1,atom1) + fx
260			f_Row(2,atom1) = f_Row(2,atom1) + fy
261			f_Row(3,atom1) = f_Row(3,atom1) + fz
262
263			f_Col(1,atom2) = f_Col(1,atom2) - fx
264			f_Col(2,atom2) = f_Col(2,atom2) - fy
265			f_Col(3,atom2) = f_Col(3,atom2) - fz
266
267	mmeineke	377	#else
268	gezelter	1160	if (do_pot) pot = pot + sw*pot_temp
269	mmeineke	377
270	gezelter	1160	f(1,atom1) = f(1,atom1) + fx
271			f(2,atom1) = f(2,atom1) + fy
272			f(3,atom1) = f(3,atom1) + fz
273
274			f(1,atom2) = f(1,atom2) - fx
275			f(2,atom2) = f(2,atom2) - fy
276			f(3,atom2) = f(3,atom2) - fz
277	mmeineke	377	#endif
278	gezelter	1192
279	gezelter	490	#ifdef IS_MPI
280	tim	1198	id1 = AtomRowToGlobal(atom1)
281			id2 = AtomColToGlobal(atom2)
282	gezelter	490	#else
283	gezelter	1192	id1 = atom1
284			id2 = atom2
285	gezelter	490	#endif
286	gezelter	1192
287			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
288	gezelter	1160
289	gezelter	1192	fpair(1) = fpair(1) + fx
290			fpair(2) = fpair(2) + fy
291			fpair(3) = fpair(3) + fz
292
293	mmeineke	377	endif
294	gezelter	1192
295	mmeineke	377	return
296	gezelter	1160
297	mmeineke	377	end subroutine do_lj_pair
298	gezelter	1160
299
300	mmeineke	377	!! Calculates the mixing for sigma or epslon
301	gezelter	1160
302	mmeineke	377	function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
303			character(len=*) :: thisParam
304			real(kind = dp) :: param1
305			real(kind = dp) :: param2
306			real(kind = dp ) :: myMixParam
307
308			integer, optional :: status
309
310			myMixParam = 0.0_dp
311
312			if (present(status)) status = 0
313			select case (LJ_Mixing_Policy)
314	gezelter	829	case (1)
315	mmeineke	377	select case (thisParam)
316			case ("sigma")
317			myMixParam = 0.5_dp * (param1 + param2)
318			case ("epsilon")
319			myMixParam = sqrt(param1 * param2)
320			case default
321			status = -1
322			end select
323			case default
324			status = -1
325			end select
326			end function calcLJMix
327
328			end module lj