UseTheForce/DarkSide/LJ.F90

!! Calculates Long Range forces Lennard-Jones interactions.
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @version $Id: LJ.F90,v 1.3 2004-10-21 20:15:25 gezelter Exp $, $Date: 2004-10-21 20:15:25 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $

module lj
  use atype_module
  use switcheroo
  use vector_class
  use simulation
  use status
#ifdef IS_MPI
  use mpiSimulation
#endif
  use force_globals

  implicit none
  PRIVATE
  
  integer, parameter :: DP = selected_real_kind(15)
  
  type, private :: LjType
     integer :: ident
     real(kind=dp) :: sigma
     real(kind=dp) :: epsilon
  end type LjType
  
  type(LjType), dimension(:), allocatable :: ParameterMap
  
  logical, save :: haveMixingMap = .false.
  
  type :: MixParameters
     real(kind=DP) :: sigma
     real(kind=DP) :: epsilon
     real(kind=dp)  :: sigma6
     real(kind=dp)  :: tp6
     real(kind=dp)  :: tp12
     real(kind=dp)  :: delta 
  end type MixParameters
  
  type(MixParameters), dimension(:,:), allocatable :: MixingMap
  
  real(kind=DP), save :: LJ_rcut
  logical, save :: have_rcut = .false.
  logical, save :: LJ_do_shift = .false.
  logical, save :: useGeometricDistanceMixing = .false.
   
  !! Public methods and data
  
  public :: setCutoffLJ
  public :: useGeometricMixing
  public :: do_lj_pair
  public :: newLJtype  
  
contains

  subroutine newLJtype(ident, sigma, epsilon, status)
    integer,intent(in) :: ident
    real(kind=dp),intent(in) :: sigma
    real(kind=dp),intent(in) :: epsilon
    integer,intent(out) :: status
    integer :: nAtypes

    status = 0
    
    !! Be simple-minded and assume that we need a ParameterMap that
    !! is the same size as the total number of atom types

    if (.not.allocated(ParameterMap)) then
       
       nAtypes = getSize(atypes)
    
       if (nAtypes == 0) then
          status = -1
          return
       end if
       
       if (.not. allocated(ParameterMap)) then
          allocate(ParameterMap(nAtypes))
       endif
       
    end if

    if (ident .gt. size(ParameterMap)) then
       status = -1
       return
    endif
    
    ! set the values for ParameterMap for this atom type:

    ParameterMap(ident)%ident = ident
    ParameterMap(ident)%epsilon = epsilon
    ParameterMap(ident)%sigma = sigma
    
  end subroutine newLJtype
    
  subroutine setCutoffLJ(rcut, do_shift, status)
    logical, intent(in):: do_shift
    integer :: status, myStatus
    real(kind=dp) :: rcut

#define __FORTRAN90
#include "UseTheForce/fSwitchingFunction.h"

    status = 0

    LJ_rcut = rcut
    LJ_do_shift = do_shift
    call set_switch(LJ_SWITCH, rcut, rcut)
    have_rcut = .true.
    
    return
  end subroutine setCutoffLJ

  subroutine useGeometricMixing() 
    useGeometricDistanceMixing = .true.
    haveMixingMap = .false.
    return
  end subroutine useGeometricMixing
  
  subroutine createMixingMap(status)
    integer :: nAtypes
    integer :: status
    integer :: i
    integer :: j
    real ( kind = dp ) :: Sigma_i, Sigma_j
    real ( kind = dp ) :: Epsilon_i, Epsilon_j
    real ( kind = dp ) :: rcut6

    status = 0
    
    nAtypes = size(ParameterMap)
    
    if (nAtypes == 0) then
       status = -1
       return
    end if

    if (.not.have_rcut) then
       status = -1
       return
    endif
    
    if (.not. allocated(MixingMap)) then
       allocate(MixingMap(nAtypes, nAtypes))
    endif
    
    rcut6 = LJ_rcut**6
    
    ! This loops through all atypes, even those that don't support LJ forces.
    do i = 1, nAtypes
       
       Epsilon_i = ParameterMap(i)%epsilon
       Sigma_i = ParameterMap(i)%sigma
       
       ! do self mixing rule
       MixingMap(i,i)%sigma   = Sigma_i          
       MixingMap(i,i)%sigma6  = Sigma_i ** 6          
       MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
       MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
       MixingMap(i,i)%epsilon = Epsilon_i          
       MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
            (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
       
       do j = i + 1, nAtypes
          
          Epsilon_j = ParameterMap(j)%epsilon
          Sigma_j = ParameterMap(j)%sigma
          
          ! only the distance parameter uses different mixing policies
          if (useGeometricDistanceMixing) then
             ! only for OPLS as far as we can tell
             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
          else
             ! everyone else
             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
          endif
          
          ! energy parameter is always geometric mean:
          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
                    
          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
          
          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
          
          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
          
       end do
    end do
    
  end subroutine createMixingMap
        
  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, localError

    if (.not.haveMixingMap) then
       localError = 0
       call createMixingMap(localError)
       if ( localError .ne. 0 ) then
          call handleError("LJ", "MixingMap creation failed!")
          return
       end if
    endif

    ! Look up the correct parameters in the mixing matrix
#ifdef IS_MPI
    sigma6   = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
    epsilon  = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
    delta    = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
#else
    sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
    epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
    delta    = MixingMap(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
    
end module lj

subroutine newLJtype(ident, sigma, epsilon, status)
  use lj, ONLY : module_newLJtype => newLJtype
  integer, parameter :: DP = selected_real_kind(15)
  integer,intent(inout) :: ident
  real(kind=dp),intent(inout) :: sigma
  real(kind=dp),intent(inout) :: epsilon
  integer,intent(inout) :: status
  
  call module_newLJtype(ident, sigma, epsilon, status)
  
end subroutine newLJtype

subroutine useGeometricMixing()
  use lj, ONLY: module_useGeometricMixing => useGeometricMixing
  
  call module_useGeometricMixing()
  return
end subroutine useGeometricMixing
Revision:	135
Committed:	Thu Oct 21 20:15:31 2004 UTC (20 years, 9 months ago) by gezelter
File size:	9093 byte(s)
Log Message:	Breaky Breaky. Fixey Fixey.
#	User	Rev	Content
1	gezelter	115	!! Calculates Long Range forces Lennard-Jones interactions.
2			!! @author Charles F. Vardeman II
3			!! @author Matthew Meineke
4	gezelter	135	!! @version $Id: LJ.F90,v 1.3 2004-10-21 20:15:25 gezelter Exp $, $Date: 2004-10-21 20:15:25 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $
5	gezelter	115
6			module lj
7			use atype_module
8			use switcheroo
9			use vector_class
10			use simulation
11	gezelter	135	use status
12	gezelter	115	#ifdef IS_MPI
13			use mpiSimulation
14			#endif
15			use force_globals
16
17			implicit none
18			PRIVATE
19	chuckv	131
20			integer, parameter :: DP = selected_real_kind(15)
21	gezelter	115
22	gezelter	135	type, private :: LjType
23			integer :: ident
24			real(kind=dp) :: sigma
25			real(kind=dp) :: epsilon
26			end type LjType
27	gezelter	115
28	gezelter	135	type(LjType), dimension(:), allocatable :: ParameterMap
29	gezelter	115
30	gezelter	135	logical, save :: haveMixingMap = .false.
31	gezelter	115
32	gezelter	135	type :: MixParameters
33			real(kind=DP) :: sigma
34			real(kind=DP) :: epsilon
35			real(kind=dp) :: sigma6
36			real(kind=dp) :: tp6
37			real(kind=dp) :: tp12
38			real(kind=dp) :: delta
39			end type MixParameters
40	chuckv	131
41	gezelter	135	type(MixParameters), dimension(:,:), allocatable :: MixingMap
42	chuckv	131
43	gezelter	135	real(kind=DP), save :: LJ_rcut
44			logical, save :: have_rcut = .false.
45			logical, save :: LJ_do_shift = .false.
46			logical, save :: useGeometricDistanceMixing = .false.
47
48			!! Public methods and data
49	chuckv	131
50	gezelter	135	public :: setCutoffLJ
51			public :: useGeometricMixing
52			public :: do_lj_pair
53			public :: newLJtype
54	chuckv	131
55	gezelter	115	contains
56
57	gezelter	135	subroutine newLJtype(ident, sigma, epsilon, status)
58	chuckv	131	integer,intent(in) :: ident
59	gezelter	135	real(kind=dp),intent(in) :: sigma
60			real(kind=dp),intent(in) :: epsilon
61	chuckv	131	integer,intent(out) :: status
62			integer :: nAtypes
63	gezelter	135
64	chuckv	131	status = 0
65
66	gezelter	135	!! Be simple-minded and assume that we need a ParameterMap that
67			!! is the same size as the total number of atom types
68	chuckv	131
69	gezelter	135	if (.not.allocated(ParameterMap)) then
70
71			nAtypes = getSize(atypes)
72
73	chuckv	131	if (nAtypes == 0) then
74	gezelter	135	status = -1
75			return
76	chuckv	131	end if
77	gezelter	135
78			if (.not. allocated(ParameterMap)) then
79			allocate(ParameterMap(nAtypes))
80			endif
81
82	chuckv	131	end if
83
84	gezelter	135	if (ident .gt. size(ParameterMap)) then
85			status = -1
86			return
87			endif
88	chuckv	131
89	gezelter	135	! set the values for ParameterMap for this atom type:
90
91			ParameterMap(ident)%ident = ident
92			ParameterMap(ident)%epsilon = epsilon
93			ParameterMap(ident)%sigma = sigma
94	chuckv	131
95			end subroutine newLJtype
96	gezelter	115
97			subroutine setCutoffLJ(rcut, do_shift, status)
98			logical, intent(in):: do_shift
99			integer :: status, myStatus
100			real(kind=dp) :: rcut
101
102			#define __FORTRAN90
103			#include "UseTheForce/fSwitchingFunction.h"
104
105			status = 0
106
107			LJ_rcut = rcut
108			LJ_do_shift = do_shift
109			call set_switch(LJ_SWITCH, rcut, rcut)
110	gezelter	135	have_rcut = .true.
111	gezelter	115
112			return
113			end subroutine setCutoffLJ
114	gezelter	135
115			subroutine useGeometricMixing()
116			useGeometricDistanceMixing = .true.
117			haveMixingMap = .false.
118			return
119			end subroutine useGeometricMixing
120	gezelter	115
121	gezelter	135	subroutine createMixingMap(status)
122	gezelter	115	integer :: nAtypes
123			integer :: status
124			integer :: i
125			integer :: j
126	gezelter	135	real ( kind = dp ) :: Sigma_i, Sigma_j
127			real ( kind = dp ) :: Epsilon_i, Epsilon_j
128	gezelter	115	real ( kind = dp ) :: rcut6
129	gezelter	135
130	gezelter	115	status = 0
131
132	gezelter	135	nAtypes = size(ParameterMap)
133
134	gezelter	115	if (nAtypes == 0) then
135			status = -1
136			return
137			end if
138	gezelter	135
139			if (.not.have_rcut) then
140			status = -1
141			return
142	gezelter	115	endif
143	gezelter	135
144			if (.not. allocated(MixingMap)) then
145			allocate(MixingMap(nAtypes, nAtypes))
146			endif
147
148	gezelter	115	rcut6 = LJ_rcut**6
149	gezelter	135
150			! This loops through all atypes, even those that don't support LJ forces.
151	gezelter	115	do i = 1, nAtypes
152	gezelter	135
153			Epsilon_i = ParameterMap(i)%epsilon
154			Sigma_i = ParameterMap(i)%sigma
155
156			! do self mixing rule
157			MixingMap(i,i)%sigma = Sigma_i
158			MixingMap(i,i)%sigma6 = Sigma_i ** 6
159			MixingMap(i,i)%tp6 = (MixingMap(i,i)%sigma6)/rcut6
160			MixingMap(i,i)%tp12 = (MixingMap(i,i)%tp6) ** 2
161			MixingMap(i,i)%epsilon = Epsilon_i
162			MixingMap(i,i)%delta = -4.0_DP * MixingMap(i,i)%epsilon * &
163			(MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
164
165			do j = i + 1, nAtypes
166	chuckv	131
167	gezelter	135	Epsilon_j = ParameterMap(j)%epsilon
168			Sigma_j = ParameterMap(j)%sigma
169	gezelter	115
170	gezelter	135	! only the distance parameter uses different mixing policies
171			if (useGeometricDistanceMixing) then
172			! only for OPLS as far as we can tell
173			MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
174			else
175			! everyone else
176			MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
177			endif
178	gezelter	115
179	gezelter	135	! energy parameter is always geometric mean:
180			MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
181
182			MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
183			MixingMap(i,j)%tp6 = MixingMap(i,j)%sigma6/rcut6
184			MixingMap(i,j)%tp12 = (MixingMap(i,j)%tp6) ** 2
185	gezelter	115
186	gezelter	135	MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
187			(MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
188	gezelter	115
189	gezelter	135	MixingMap(j,i)%sigma = MixingMap(i,j)%sigma
190			MixingMap(j,i)%sigma6 = MixingMap(i,j)%sigma6
191			MixingMap(j,i)%tp6 = MixingMap(i,j)%tp6
192			MixingMap(j,i)%tp12 = MixingMap(i,j)%tp12
193			MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
194			MixingMap(j,i)%delta = MixingMap(i,j)%delta
195	gezelter	115
196	gezelter	135	end do
197	gezelter	115	end do
198
199	gezelter	135	end subroutine createMixingMap
200
201	gezelter	115	subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
202			pot, f, do_pot)
203
204			integer, intent(in) :: atom1, atom2
205			real( kind = dp ), intent(in) :: rij, r2
206			real( kind = dp ) :: pot, sw, vpair
207			real( kind = dp ), dimension(3,nLocal) :: f
208			real( kind = dp ), intent(in), dimension(3) :: d
209			real( kind = dp ), intent(inout), dimension(3) :: fpair
210			logical, intent(in) :: do_pot
211
212			! local Variables
213			real( kind = dp ) :: drdx, drdy, drdz
214			real( kind = dp ) :: fx, fy, fz
215			real( kind = dp ) :: pot_temp, dudr
216			real( kind = dp ) :: sigma6
217			real( kind = dp ) :: epsilon
218			real( kind = dp ) :: r6
219			real( kind = dp ) :: t6
220			real( kind = dp ) :: t12
221			real( kind = dp ) :: delta
222	gezelter	135	integer :: id1, id2, localError
223	gezelter	115
224	gezelter	135	if (.not.haveMixingMap) then
225			localError = 0
226			call createMixingMap(localError)
227			if ( localError .ne. 0 ) then
228			call handleError("LJ", "MixingMap creation failed!")
229			return
230			end if
231			endif
232
233	gezelter	115	! Look up the correct parameters in the mixing matrix
234			#ifdef IS_MPI
235	gezelter	135	sigma6 = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
236			epsilon = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
237			delta = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
238	gezelter	115	#else
239	gezelter	135	sigma6 = MixingMap(atid(atom1),atid(atom2))%sigma6
240			epsilon = MixingMap(atid(atom1),atid(atom2))%epsilon
241			delta = MixingMap(atid(atom1),atid(atom2))%delta
242	gezelter	115	#endif
243
244			r6 = r2 * r2 * r2
245
246			t6 = sigma6/ r6
247			t12 = t6 * t6
248
249			pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
250			if (LJ_do_shift) then
251			pot_temp = pot_temp + delta
252			endif
253
254			vpair = vpair + pot_temp
255
256			dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
257
258			drdx = d(1) / rij
259			drdy = d(2) / rij
260			drdz = d(3) / rij
261
262			fx = dudr * drdx
263			fy = dudr * drdy
264			fz = dudr * drdz
265
266
267			#ifdef IS_MPI
268			if (do_pot) then
269			pot_Row(atom1) = pot_Row(atom1) + swpot_temp0.5
270			pot_Col(atom2) = pot_Col(atom2) + swpot_temp0.5
271			endif
272
273			f_Row(1,atom1) = f_Row(1,atom1) + fx
274			f_Row(2,atom1) = f_Row(2,atom1) + fy
275			f_Row(3,atom1) = f_Row(3,atom1) + fz
276
277			f_Col(1,atom2) = f_Col(1,atom2) - fx
278			f_Col(2,atom2) = f_Col(2,atom2) - fy
279			f_Col(3,atom2) = f_Col(3,atom2) - fz
280
281			#else
282			if (do_pot) pot = pot + sw*pot_temp
283
284			f(1,atom1) = f(1,atom1) + fx
285			f(2,atom1) = f(2,atom1) + fy
286			f(3,atom1) = f(3,atom1) + fz
287
288			f(1,atom2) = f(1,atom2) - fx
289			f(2,atom2) = f(2,atom2) - fy
290			f(3,atom2) = f(3,atom2) - fz
291			#endif
292
293			#ifdef IS_MPI
294			id1 = AtomRowToGlobal(atom1)
295			id2 = AtomColToGlobal(atom2)
296			#else
297			id1 = atom1
298			id2 = atom2
299			#endif
300
301			if (molMembershipList(id1) .ne. molMembershipList(id2)) then
302
303			fpair(1) = fpair(1) + fx
304			fpair(2) = fpair(2) + fy
305			fpair(3) = fpair(3) + fz
306
307			endif
308
309			return
310
311			end subroutine do_lj_pair
312
313
314			!! Calculates the mixing for sigma or epslon
315
316			end module lj
317	chuckv	131
318	gezelter	135	subroutine newLJtype(ident, sigma, epsilon, status)
319			use lj, ONLY : module_newLJtype => newLJtype
320			integer, parameter :: DP = selected_real_kind(15)
321			integer,intent(inout) :: ident
322			real(kind=dp),intent(inout) :: sigma
323			real(kind=dp),intent(inout) :: epsilon
324			integer,intent(inout) :: status
325
326			call module_newLJtype(ident, sigma, epsilon, status)
327
328			end subroutine newLJtype
329	chuckv	131
330	gezelter	135	subroutine useGeometricMixing()
331			use lj, ONLY: module_useGeometricMixing => useGeometricMixing
332
333			call module_useGeometricMixing()
334			return
335			end subroutine useGeometricMixing