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 1628 by gezelter, Thu Oct 21 20:15:31 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.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  
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    
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  
55   contains
56  
57 <  subroutine newLJtype(ident,lj_sigma,lj_epsilon,status)
57 >  subroutine newLJtype(ident, sigma, epsilon, status)
58      integer,intent(in) :: ident
59 <    real(kind=dp),intent(in) :: lj_sigma
60 <    real(kind=dp),intent(in) :: lj_epsilon
59 >    real(kind=dp),intent(in) :: sigma
60 >    real(kind=dp),intent(in) :: epsilon
61      integer,intent(out) :: status
83    
84    integer,pointer                        :: Matchlist(:) => null()
85    integer :: current
62      integer :: nAtypes
63 +
64      status = 0
65      
66 <        !! Assume that atypes has already been set and get the total number of types in atypes
67 <  
91 <  
66 >    !! Be simple-minded and assume that we need a ParameterMap that
67 >    !! is the same size as the total number of atom types
68  
69 <    ! check to see if this is the first time into
70 <    if (.not.associated(ljParameterList%ljParams)) then
71 <       call getMatchingElementList(atypes, "is_lj", .true., nAtypes, MatchList)
72 <       ljParameterList%n_lj_types = nAtypes
69 >    if (.not.allocated(ParameterMap)) then
70 >      
71 >       nAtypes = getSize(atypes)
72 >    
73         if (nAtypes == 0) then
74 <         status = -1
75 <         return
74 >          status = -1
75 >          return
76         end if
77 <       allocate(ljParameterList%ljParams(nAtypes))
77 >      
78 >       if (.not. allocated(ParameterMap)) then
79 >          allocate(ParameterMap(nAtypes))
80 >       endif
81 >      
82      end if
83  
84 <    ljParameterList%currentAddition = ljParameterList%currentAddition + 1
85 <    current = ljParameterList%currentAddition
84 >    if (ident .gt. size(ParameterMap)) then
85 >       status = -1
86 >       return
87 >    endif
88      
89 <    ! set the values for ljParameterList
90 <    ljParameterList%ljParams(current)%lj_ident = ident
91 <    ljParameterList%ljParams(current)%lj_epsilon = lj_epsilon
92 <    ljParameterList%ljParams(current)%lj_sigma = lj_sigma
89 >    ! 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      
95    end subroutine newLJtype
113  
114  subroutine init_LJ_FF(mix_Policy, status)
115    integer, intent(in) :: mix_Policy
116    integer, intent(out) :: status
117    integer :: myStatus
96      
119    if (mix_Policy == LB_MIXING_RULE) then
120       LJ_Mixing_Policy = LB_MIXING_RULE
121    else
122       if (mix_Policy == EXPLICIT_MIXING_RULE) then
123          LJ_Mixing_Policy = EXPLICIT_MIXING_RULE
124       else
125          write(*,*) 'Unknown Mixing Policy!'
126          status = -1
127          return
128       endif
129    endif
130
131    havePolicy = .true.
132
133    if (haveCut) then
134       status = 0
135       call createMixingList(myStatus)
136       if (myStatus /= 0) then
137          status = -1
138          return
139       end if
140      
141       LJ_FF_initialized = .true.
142    end if
143  
144  end subroutine init_LJ_FF
145  
97    subroutine setCutoffLJ(rcut, do_shift, status)
98      logical, intent(in):: do_shift
99      integer :: status, myStatus
# Line 156 | Line 107 | contains
107      LJ_rcut = rcut
108      LJ_do_shift = do_shift
109      call set_switch(LJ_SWITCH, rcut, rcut)
110 <    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    
110 >    have_rcut = .true.
111      
112      return
113    end subroutine setCutoffLJ
114 +
115 +  subroutine useGeometricMixing()
116 +    useGeometricDistanceMixing = .true.
117 +    haveMixingMap = .false.
118 +    return
119 +  end subroutine useGeometricMixing
120    
121 <  subroutine createMixingList(status)
121 >  subroutine createMixingMap(status)
122      integer :: nAtypes
123      integer :: status
124      integer :: i
125      integer :: j
126 <    real ( kind = dp ) :: mySigma_i,mySigma_j
127 <    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
126 >    real ( kind = dp ) :: Sigma_i, Sigma_j
127 >    real ( kind = dp ) :: Epsilon_i, Epsilon_j
128      real ( kind = dp ) :: rcut6
129 <    logical :: I_isLJ, J_isLJ
129 >
130      status = 0
131      
132 <    ! we only allocate this array to the number of lj_atypes
133 <    nAtypes = size(ljParameterList%ljParams)
132 >    nAtypes = size(ParameterMap)
133 >    
134      if (nAtypes == 0) then
135         status = -1
136         return
137      end if
192        
193    if (.not. associated(ljMixed)) then
194       allocate(ljMixed(nAtypes, nAtypes))
195    endif
138  
139 <    rcut6 = LJ_rcut**6
140 <
141 < ! This loops through all atypes, even those that don't support LJ forces.
139 >    if (.not.have_rcut) then
140 >       status = -1
141 >       return
142 >    endif
143 >    
144 >    if (.not. allocated(MixingMap)) then
145 >       allocate(MixingMap(nAtypes, nAtypes))
146 >    endif
147 >    
148 >    rcut6 = LJ_rcut**6
149 >    
150 >    ! This loops through all atypes, even those that don't support LJ forces.
151      do i = 1, nAtypes
152 <
153 <          myEpsilon_i = ljParameterList%ljParams(i)%lj_epsilon
154 <          mySigma_i = ljParameterList%ljParams(i)%lj_sigma
152 >      
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            
167 <          ! do self mixing rule
168 <          ljMixed(i,i)%sigma   = mySigma_i
167 >          Epsilon_j = ParameterMap(j)%epsilon
168 >          Sigma_j = ParameterMap(j)%sigma
169            
170 <          ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6
170 >          ! 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            
179 <          ljMixed(i,i)%tp6     = (ljMixed(i,i)%sigma6)/rcut6
179 >          ! 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            
186 <          ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
186 >          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
187 >               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
188            
189 +          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            
196 <          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
196 >       end do
197      end do
198      
199 <  end subroutine createMixingList
200 <  
199 >  end subroutine createMixingMap
200 >        
201    subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
202         pot, f, do_pot)
203  
# Line 277 | Line 219 | contains
219      real( kind = dp ) :: t6
220      real( kind = dp ) :: t12
221      real( kind = dp ) :: delta
222 <    integer :: id1, id2
222 >    integer :: id1, id2, localError
223  
224 +    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      ! Look up the correct parameters in the mixing matrix
234   #ifdef IS_MPI
235 <    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
236 <    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
237 <    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
235 >    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   #else
239 <    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
240 <    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
241 <    delta    = ljMixed(atid(atom1),atid(atom2))%delta
239 >    sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
240 >    epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
241 >    delta    = MixingMap(atid(atom1),atid(atom2))%delta
242   #endif
243  
244      r6 = r2 * r2 * r2
# Line 361 | Line 312 | contains
312    
313    
314    !! 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
315      
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  
316   end module lj
317  
318 < subroutine newLJtype(ident,lj_sigma,lj_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) :: lj_sigma
323 <    real(kind=dp),intent(inout) :: lj_epsilon
324 <    integer,intent(inout) :: status
318 > 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  
330 <    call module_newLJtype(ident,lj_sigma,lj_epsilon,status)
331 <
332 < end subroutine newLJtype
333 <
330 > subroutine useGeometricMixing()
331 >  use lj, ONLY: module_useGeometricMixing => useGeometricMixing
332 >  
333 >  call module_useGeometricMixing()
334 >  return
335 > end subroutine useGeometricMixing

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines