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:	1628
Committed:	Thu Oct 21 20:15:31 2004 UTC (19 years, 10 months ago) by gezelter
File size:	9093 byte(s)
Log Message:	Breaky Breaky. Fixey Fixey.
#	User	Rev	Content
1	gezelter	1608	!! Calculates Long Range forces Lennard-Jones interactions.
2			!! @author Charles F. Vardeman II
3			!! @author Matthew Meineke
4	gezelter	1628	!! @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	1608
6			module lj
7			use atype_module
8			use switcheroo
9			use vector_class
10			use simulation
11	gezelter	1628	use status
12	gezelter	1608	#ifdef IS_MPI
13			use mpiSimulation
14			#endif
15			use force_globals
16
17			implicit none
18			PRIVATE
19	chuckv	1624
20			integer, parameter :: DP = selected_real_kind(15)
21	gezelter	1608
22	gezelter	1628	type, private :: LjType
23			integer :: ident
24			real(kind=dp) :: sigma
25			real(kind=dp) :: epsilon
26			end type LjType
27	gezelter	1608
28	gezelter	1628	type(LjType), dimension(:), allocatable :: ParameterMap
29	gezelter	1608
30	gezelter	1628	logical, save :: haveMixingMap = .false.
31	gezelter	1608
32	gezelter	1628	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	1624
41	gezelter	1628	type(MixParameters), dimension(:,:), allocatable :: MixingMap
42	chuckv	1624
43	gezelter	1628	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	1624
50	gezelter	1628	public :: setCutoffLJ
51			public :: useGeometricMixing
52			public :: do_lj_pair
53			public :: newLJtype
54	chuckv	1624
55	gezelter	1608	contains
56
57	gezelter	1628	subroutine newLJtype(ident, sigma, epsilon, status)
58	chuckv	1624	integer,intent(in) :: ident
59	gezelter	1628	real(kind=dp),intent(in) :: sigma
60			real(kind=dp),intent(in) :: epsilon
61	chuckv	1624	integer,intent(out) :: status
62			integer :: nAtypes
63	gezelter	1628
64	chuckv	1624	status = 0
65
66	gezelter	1628	!! 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	1624
69	gezelter	1628	if (.not.allocated(ParameterMap)) then
70
71			nAtypes = getSize(atypes)
72
73	chuckv	1624	if (nAtypes == 0) then
74	gezelter	1628	status = -1
75			return
76	chuckv	1624	end if
77	gezelter	1628
78			if (.not. allocated(ParameterMap)) then
79			allocate(ParameterMap(nAtypes))
80			endif
81
82	chuckv	1624	end if
83
84	gezelter	1628	if (ident .gt. size(ParameterMap)) then
85			status = -1
86			return
87			endif
88	chuckv	1624
89	gezelter	1628	! 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	1624
95			end subroutine newLJtype
96	gezelter	1608
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	1628	have_rcut = .true.
111	gezelter	1608
112			return
113			end subroutine setCutoffLJ
114	gezelter	1628
115			subroutine useGeometricMixing()
116			useGeometricDistanceMixing = .true.
117			haveMixingMap = .false.
118			return
119			end subroutine useGeometricMixing
120	gezelter	1608
121	gezelter	1628	subroutine createMixingMap(status)
122	gezelter	1608	integer :: nAtypes
123			integer :: status
124			integer :: i
125			integer :: j
126	gezelter	1628	real ( kind = dp ) :: Sigma_i, Sigma_j
127			real ( kind = dp ) :: Epsilon_i, Epsilon_j
128	gezelter	1608	real ( kind = dp ) :: rcut6
129	gezelter	1628
130	gezelter	1608	status = 0
131
132	gezelter	1628	nAtypes = size(ParameterMap)
133
134	gezelter	1608	if (nAtypes == 0) then
135			status = -1
136			return
137			end if
138	gezelter	1628
139			if (.not.have_rcut) then
140			status = -1
141			return
142	gezelter	1608	endif
143	gezelter	1628
144			if (.not. allocated(MixingMap)) then
145			allocate(MixingMap(nAtypes, nAtypes))
146			endif
147
148	gezelter	1608	rcut6 = LJ_rcut**6
149	gezelter	1628
150			! This loops through all atypes, even those that don't support LJ forces.
151	gezelter	1608	do i = 1, nAtypes
152	gezelter	1628
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	1624
167	gezelter	1628	Epsilon_j = ParameterMap(j)%epsilon
168			Sigma_j = ParameterMap(j)%sigma
169	gezelter	1608
170	gezelter	1628	! 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	1608
179	gezelter	1628	! 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	1608
186	gezelter	1628	MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
187			(MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
188	gezelter	1608
189	gezelter	1628	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	1608
196	gezelter	1628	end do
197	gezelter	1608	end do
198
199	gezelter	1628	end subroutine createMixingMap
200
201	gezelter	1608	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	1628	integer :: id1, id2, localError
223	gezelter	1608
224	gezelter	1628	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	1608	! Look up the correct parameters in the mixing matrix
234			#ifdef IS_MPI
235	gezelter	1628	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	1608	#else
239	gezelter	1628	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	1608	#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	1624
318	gezelter	1628	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	1624
330	gezelter	1628	subroutine useGeometricMixing()
331			use lj, ONLY: module_useGeometricMixing => useGeometricMixing
332
333			call module_useGeometricMixing()
334			return
335			end subroutine useGeometricMixing