ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/LJ.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/LJ.F90 (file contents):
Revision 1624 by chuckv, Thu Oct 21 15:25:30 2004 UTC vs.
Revision 1633 by gezelter, Fri Oct 22 20:22:48 2004 UTC

# Line 1 | Line 1
1   !! Calculates Long Range forces Lennard-Jones interactions.
2 !! Corresponds to the force field defined in lj_FF.cpp
2   !! @author Charles F. Vardeman II
3   !! @author Matthew Meineke
4 < !! @version $Id: LJ.F90,v 1.2 2004-10-21 15:25:30 chuckv Exp $, $Date: 2004-10-21 15:25:30 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
4 > !! @version $Id: LJ.F90,v 1.4 2004-10-22 20:22:47 gezelter Exp $, $Date: 2004-10-22 20:22:47 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
5  
6   module lj
7    use atype_module
8    use switcheroo
9    use vector_class
10    use simulation
11 +  use status
12   #ifdef IS_MPI
13    use mpiSimulation
14   #endif
# Line 18 | Line 18 | module lj
18    PRIVATE
19    
20    integer, parameter :: DP = selected_real_kind(15)
21
22 #define __FORTRAN90
23 #include "UseTheForce/fForceField.h"
24
25  integer, save :: LJ_Mixing_Policy
26  real(kind=DP), save :: LJ_rcut
27  logical, save :: havePolicy = .false.
28  logical, save :: haveCut = .false.
29  logical, save :: LJ_do_shift = .false.
21    
22 <  !! Logical has lj force field module been initialized?
22 >  type, private :: LjType
23 >     integer :: ident
24 >     real(kind=dp) :: sigma
25 >     real(kind=dp) :: epsilon
26 >  end type LjType
27    
28 <  logical, save :: LJ_FF_initialized = .false.
28 >  type(LjType), dimension(:), allocatable :: ParameterMap
29    
30 <  !! Public methods and data
36 <  public :: init_LJ_FF
37 <  public :: setCutoffLJ
38 <  public :: do_lj_pair
39 <  public :: newLJtype
30 >  logical, save :: haveMixingMap = .false.
31    
32 <  !! structure for lj type parameters
33 <  type, private :: ljType
34 <    integer :: lj_ident
35 <    real(kind=dp) :: lj_sigma
36 <    real(kind=dp) :: lj_epsilon
37 <  end type ljType
32 >  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    
41 <  !! List of lj type parameters
49 <  type, private :: ljTypeList
50 <    integer  :: n_lj_types = 0
51 <    integer  :: currentAddition = 0
52 <    type(ljType), pointer :: ljParams(:) => null()
53 <  end type ljTypeList
41 >  type(MixParameters), dimension(:,:), allocatable :: MixingMap
42    
43 <  !! The list of lj Parameters
44 <  type (ljTypeList), save :: ljParameterList
43 >  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    
50 +  public :: setCutoffLJ
51 +  public :: useGeometricMixing
52 +  public :: do_lj_pair
53 +  public :: newLJtype  
54 +  public :: getSigma
55 +  public :: getEpsilon
56    
59  type :: lj_mixed_params
60     !! Lennard-Jones epsilon
61     real ( kind = dp )  :: epsilon = 0.0_dp
62     !! Lennard-Jones Sigma
63     real ( kind = dp )  :: sigma = 0.0_dp
64     !! Lennard-Jones Sigma to sixth
65     real ( kind = dp )  :: sigma6 = 0.0_dp
66     !!
67     real ( kind = dp )  :: tp6
68     real ( kind = dp )  :: tp12
69     real ( kind = dp )  :: delta  = 0.0_dp
70  end type lj_mixed_params
71  
72  type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
73  
74  
75  
57   contains
58  
59 <  subroutine newLJtype(ident,lj_sigma,lj_epsilon,status)
59 >  subroutine newLJtype(ident, sigma, epsilon, status)
60      integer,intent(in) :: ident
61 <    real(kind=dp),intent(in) :: lj_sigma
62 <    real(kind=dp),intent(in) :: lj_epsilon
61 >    real(kind=dp),intent(in) :: sigma
62 >    real(kind=dp),intent(in) :: epsilon
63      integer,intent(out) :: status
83    
84    integer,pointer                        :: Matchlist(:) => null()
85    integer :: current
64      integer :: nAtypes
65 +
66      status = 0
67      
68 <        !! Assume that atypes has already been set and get the total number of types in atypes
69 <  
91 <  
68 >    !! Be simple-minded and assume that we need a ParameterMap that
69 >    !! is the same size as the total number of atom types
70  
71 <    ! check to see if this is the first time into
72 <    if (.not.associated(ljParameterList%ljParams)) then
73 <       call getMatchingElementList(atypes, "is_lj", .true., nAtypes, MatchList)
74 <       ljParameterList%n_lj_types = nAtypes
71 >    if (.not.allocated(ParameterMap)) then
72 >      
73 >       nAtypes = getSize(atypes)
74 >    
75         if (nAtypes == 0) then
76 <         status = -1
77 <         return
76 >          status = -1
77 >          return
78         end if
79 <       allocate(ljParameterList%ljParams(nAtypes))
79 >      
80 >       if (.not. allocated(ParameterMap)) then
81 >          allocate(ParameterMap(nAtypes))
82 >       endif
83 >      
84      end if
85  
86 <    ljParameterList%currentAddition = ljParameterList%currentAddition + 1
87 <    current = ljParameterList%currentAddition
86 >    if (ident .gt. size(ParameterMap)) then
87 >       status = -1
88 >       return
89 >    endif
90      
91 <    ! set the values for ljParameterList
92 <    ljParameterList%ljParams(current)%lj_ident = ident
93 <    ljParameterList%ljParams(current)%lj_epsilon = lj_epsilon
94 <    ljParameterList%ljParams(current)%lj_sigma = lj_sigma
91 >    ! set the values for ParameterMap for this atom type:
92 >
93 >    ParameterMap(ident)%ident = ident
94 >    ParameterMap(ident)%epsilon = epsilon
95 >    ParameterMap(ident)%sigma = sigma
96      
97    end subroutine newLJtype
98 <  
99 <  subroutine init_LJ_FF(mix_Policy, status)
100 <    integer, intent(in) :: mix_Policy
101 <    integer, intent(out) :: status
102 <    integer :: myStatus
98 >
99 >  function getSigma(atid) result (s)
100 >    integer, intent(in) :: atid
101 >    integer :: localError
102 >    real(kind=dp) :: s
103      
104 <    if (mix_Policy == LB_MIXING_RULE) then
105 <       LJ_Mixing_Policy = LB_MIXING_RULE
106 <    else
107 <       if (mix_Policy == EXPLICIT_MIXING_RULE) then
108 <          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
109 <       else
110 <          write(*,*) 'Unknown Mixing Policy!'
126 <          status = -1
127 <          return
128 <       endif
129 <    endif
104 >    if (.not.allocated(ParameterMap)) then
105 >       call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
106 >       return
107 >    end if
108 >    
109 >    s = ParameterMap(atid)%sigma
110 >  end function getSigma
111  
112 <    havePolicy = .true.
113 <
114 <    if (haveCut) then
115 <       status = 0
116 <       call createMixingList(myStatus)
117 <       if (myStatus /= 0) then
118 <          status = -1
119 <          return
139 <       end if
140 <      
141 <       LJ_FF_initialized = .true.
112 >  function getEpsilon(atid) result (e)
113 >    integer, intent(in) :: atid
114 >    integer :: localError
115 >    real(kind=dp) :: e
116 >    
117 >    if (.not.allocated(ParameterMap)) then
118 >       call handleError("dipole-dipole", "no ParameterMap was present before first call of getEpsilon!")
119 >       return
120      end if
121 <  
122 <  end subroutine init_LJ_FF
123 <  
121 >    
122 >    e = ParameterMap(atid)%epsilon
123 >  end function getEpsilon
124 >
125 >
126    subroutine setCutoffLJ(rcut, do_shift, status)
127      logical, intent(in):: do_shift
128      integer :: status, myStatus
# Line 156 | Line 136 | contains
136      LJ_rcut = rcut
137      LJ_do_shift = do_shift
138      call set_switch(LJ_SWITCH, rcut, rcut)
139 <    haveCut = .true.
160 <
161 <    if (havePolicy) then
162 <       status = 0
163 <       call createMixingList(myStatus)
164 <       if (myStatus /= 0) then
165 <          status = -1
166 <          return
167 <       end if
168 <      
169 <       LJ_FF_initialized = .true.
170 <    end if    
139 >    have_rcut = .true.
140      
141      return
142    end subroutine setCutoffLJ
143 +
144 +  subroutine useGeometricMixing()
145 +    useGeometricDistanceMixing = .true.
146 +    haveMixingMap = .false.
147 +    return
148 +  end subroutine useGeometricMixing
149    
150 <  subroutine createMixingList(status)
150 >  subroutine createMixingMap(status)
151      integer :: nAtypes
152      integer :: status
153      integer :: i
154      integer :: j
155 <    real ( kind = dp ) :: mySigma_i,mySigma_j
156 <    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
155 >    real ( kind = dp ) :: Sigma_i, Sigma_j
156 >    real ( kind = dp ) :: Epsilon_i, Epsilon_j
157      real ( kind = dp ) :: rcut6
158 <    logical :: I_isLJ, J_isLJ
158 >
159      status = 0
160      
161 <    ! we only allocate this array to the number of lj_atypes
162 <    nAtypes = size(ljParameterList%ljParams)
161 >    nAtypes = size(ParameterMap)
162 >    
163      if (nAtypes == 0) then
164         status = -1
165         return
166      end if
192        
193    if (.not. associated(ljMixed)) then
194       allocate(ljMixed(nAtypes, nAtypes))
195    endif
167  
168 +    if (.not.have_rcut) then
169 +       status = -1
170 +       return
171 +    endif
172 +    
173 +    if (.not. allocated(MixingMap)) then
174 +       allocate(MixingMap(nAtypes, nAtypes))
175 +    endif
176 +    
177      rcut6 = LJ_rcut**6
178 <
179 < ! This loops through all atypes, even those that don't support LJ forces.
178 >    
179 >    ! This loops through all atypes, even those that don't support LJ forces.
180      do i = 1, nAtypes
181 <
182 <          myEpsilon_i = ljParameterList%ljParams(i)%lj_epsilon
183 <          mySigma_i = ljParameterList%ljParams(i)%lj_sigma
181 >      
182 >       Epsilon_i = ParameterMap(i)%epsilon
183 >       Sigma_i = ParameterMap(i)%sigma
184 >      
185 >       ! do self mixing rule
186 >       MixingMap(i,i)%sigma   = Sigma_i          
187 >       MixingMap(i,i)%sigma6  = Sigma_i ** 6          
188 >       MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
189 >       MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
190 >       MixingMap(i,i)%epsilon = Epsilon_i          
191 >       MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
192 >            (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
193 >      
194 >       do j = i + 1, nAtypes
195            
196 <          ! do self mixing rule
197 <          ljMixed(i,i)%sigma   = mySigma_i
196 >          Epsilon_j = ParameterMap(j)%epsilon
197 >          Sigma_j = ParameterMap(j)%sigma
198            
199 <          ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
199 >          ! only the distance parameter uses different mixing policies
200 >          if (useGeometricDistanceMixing) then
201 >             ! only for OPLS as far as we can tell
202 >             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
203 >          else
204 >             ! everyone else
205 >             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
206 >          endif
207            
208 <          ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
208 >          ! energy parameter is always geometric mean:
209 >          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
210 >                    
211 >          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
212 >          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
213 >          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
214            
215 <          ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
215 >          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
216 >               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
217            
218 +          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
219 +          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
220 +          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
221 +          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
222 +          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
223 +          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
224            
225 <          ljMixed(i,i)%epsilon = myEpsilon_i
216 <          
217 <          ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
218 <            (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
219 <          
220 <          do j = i + 1, nAtypes
221 <
222 <                myEpsilon_j = ljParameterList%ljParams(j)%lj_epsilon
223 <                mySigma_j = ljParameterList%ljParams(j)%lj_sigma
224 <
225 <                          
226 <                ljMixed(i,j)%sigma  =  &
227 <                     calcLJMix("sigma",mySigma_i, &
228 <                     mySigma_j)
229 <                
230 <                ljMixed(i,j)%sigma6 = &
231 <                     (ljMixed(i,j)%sigma)**6
232 <                
233 <                
234 <                ljMixed(i,j)%tp6     = ljMixed(i,j)%sigma6/rcut6
235 <                
236 <                ljMixed(i,j)%tp12    = (ljMixed(i,j)%tp6) ** 2
237 <                
238 <                
239 <                ljMixed(i,j)%epsilon = &
240 <                     calcLJMix("epsilon",myEpsilon_i, &
241 <                     myEpsilon_j)
242 <                
243 <                ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
244 <                     (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
245 <                
246 <                
247 <                ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
248 <                ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
249 <                ljMixed(j,i)%tp6     = ljMixed(i,j)%tp6
250 <                ljMixed(j,i)%tp12    = ljMixed(i,j)%tp12
251 <                ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
252 <                ljMixed(j,i)%delta   = ljMixed(i,j)%delta
253 <          
254 <          end do
225 >       end do
226      end do
227      
228 <  end subroutine createMixingList
229 <  
228 >  end subroutine createMixingMap
229 >        
230    subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
231         pot, f, do_pot)
232  
# Line 277 | Line 248 | contains
248      real( kind = dp ) :: t6
249      real( kind = dp ) :: t12
250      real( kind = dp ) :: delta
251 <    integer :: id1, id2
251 >    integer :: id1, id2, localError
252  
253 +    if (.not.haveMixingMap) then
254 +       localError = 0
255 +       call createMixingMap(localError)
256 +       if ( localError .ne. 0 ) then
257 +          call handleError("LJ", "MixingMap creation failed!")
258 +          return
259 +       end if
260 +    endif
261 +
262      ! Look up the correct parameters in the mixing matrix
263   #ifdef IS_MPI
264 <    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
265 <    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
266 <    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
264 >    sigma6   = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
265 >    epsilon  = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
266 >    delta    = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
267   #else
268 <    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
269 <    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
270 <    delta    = ljMixed(atid(atom1),atid(atom2))%delta
268 >    sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
269 >    epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
270 >    delta    = MixingMap(atid(atom1),atid(atom2))%delta
271   #endif
272  
273      r6 = r2 * r2 * r2
# Line 361 | Line 341 | contains
341    
342    
343    !! Calculates the mixing for sigma or epslon
364  
365  function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
366    character(len=*) :: thisParam
367    real(kind = dp)  :: param1
368    real(kind = dp)  :: param2
369    real(kind = dp ) :: myMixParam
370
371    integer, optional :: status  
372
373    myMixParam = 0.0_dp
344      
375    if (present(status)) status = 0
376    select case (LJ_Mixing_Policy)
377    case (1)
378       select case (thisParam)
379       case ("sigma")
380          myMixParam = 0.5_dp * (param1 + param2)
381       case ("epsilon")
382          myMixParam = sqrt(param1 * param2)
383       case default
384          status = -1
385       end select
386    case default
387       status = -1
388    end select
389  end function calcLJMix
390  
345   end module lj
346  
347 < subroutine newLJtype(ident,lj_sigma,lj_epsilon,status)
348 <    use lj, ONLY : module_newLJtype => newLJtype
349 <    integer, parameter :: DP = selected_real_kind(15)
350 <    integer,intent(inout) :: ident
351 <    real(kind=dp),intent(inout) :: lj_sigma
352 <    real(kind=dp),intent(inout) :: lj_epsilon
353 <    integer,intent(inout) :: status
347 > subroutine newLJtype(ident, sigma, epsilon, status)
348 >  use lj, ONLY : module_newLJtype => newLJtype
349 >  integer, parameter :: DP = selected_real_kind(15)
350 >  integer,intent(inout) :: ident
351 >  real(kind=dp),intent(inout) :: sigma
352 >  real(kind=dp),intent(inout) :: epsilon
353 >  integer,intent(inout) :: status
354 >  
355 >  call module_newLJtype(ident, sigma, epsilon, status)
356 >  
357 > end subroutine newLJtype
358  
359 <    call module_newLJtype(ident,lj_sigma,lj_epsilon,status)
360 <
361 < end subroutine newLJtype
362 <
359 > subroutine useGeometricMixing()
360 >  use lj, ONLY: module_useGeometricMixing => useGeometricMixing
361 >  
362 >  call module_useGeometricMixing()
363 >  return
364 > end subroutine useGeometricMixing

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines