ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/calc_LJ_FF.F90
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/calc_LJ_FF.F90 (file contents):
Revision 306 by chuckv, Mon Mar 10 19:26:45 2003 UTC vs.
Revision 329 by gezelter, Wed Mar 12 22:27:59 2003 UTC

# Line 2 | Line 2
2   !! Corresponds to the force field defined in lj_FF.cpp
3   !! @author Charles F. Vardeman II
4   !! @author Matthew Meineke
5 < !! @version $Id: calc_LJ_FF.F90,v 1.2 2003-03-10 19:26:45 chuckv Exp $, $Date: 2003-03-10 19:26:45 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
5 > !! @version $Id: calc_LJ_FF.F90,v 1.9 2003-03-12 22:27:59 gezelter Exp $, $Date: 2003-03-12 22:27:59 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $
6  
7 <
8 <
9 < module lj_ff
7 > module lj
8    use simulation
9    use definitions
10 <  use generic_lists
11 <
10 >  use atype_module
11 >  use vector_class
12   #ifdef IS_MPI
13    use mpiSimulation
14   #endif
15    implicit none
16    PRIVATE
17  
18 < !! Logical has lj force field module been initialized?
19 <  logical, save :: isljFFinit = .false.
20 <
21 <
22 < !! Public methods and data
23 <  public :: getLjPot
24 <  public :: init_ljFF
25 <
18 >  integer, save :: LJ_Mixing_Policy
19 >  integer, parameter :: LB_Mixing_Rule = 1
20 >  integer, parameter :: Explicit_Mixing_Rule = 2
21 >  
22 >  !! Logical has lj force field module been initialized?
23 >  
24 >  logical, save :: LJ_FF_initialized = .false.
25 >  
26 >  !! Public methods and data
27 >  public :: init_LJ_FF
28 >  public :: do_lj_pair
29 >  
30    type :: lj_mixed_params
31 <     !! Mass of Particle
30 <     real ( kind = dp )  :: mass = 0.0_dp
31 < !! Lennard-Jones epslon
31 >     !! Lennard-Jones epsilon
32       real ( kind = dp )  :: epsilon = 0.0_dp
33 < !! Lennard-Jones Sigma
33 >     !! Lennard-Jones Sigma
34       real ( kind = dp )  :: sigma = 0.0_dp
35 < !! Lennard-Jones Sigma Squared
36 <     real ( kind = dp )  :: sigma2 = 0.0_dp
37 < !! Lennard-Jones Sigma to sixth
35 >     !! Lennard-Jones Sigma to sixth
36       real ( kind = dp )  :: sigma6 = 0.0_dp
37 +     !!
38 +     real ( kind = dp )  :: tp6
39 +     real ( kind = dp )  :: tp12
40 +
41 +     real ( kind = dp )  :: delta  = 0.0_dp
42 +
43    end type lj_mixed_params
44    
45 <  type (lj_mixed_params), dimension(:,:) :: ljMixed
46 <
45 >  type (lj_mixed_params), dimension(:,:), pointer :: ljMixed
46 >  
47   contains
48  
49 <  subroutine init_ljFF(ListHead,status)
50 <    type (atype),pointer :: ListHead
47 <    integer :: nTypes
49 >  subroutine init_LJ_FF(mix_Policy, status)
50 >    character(len=100), intent(in) :: mix_Policy
51      integer, intent(out) :: status
49
52      integer :: myStatus
53 <
53 >    
54 >    if ((mix_Policy == "standard") .or.  (mix_Policy == "LB")) then
55 >       LJ_Mixing_Policy = LB_Mixing_Rule
56 >    else
57 >       if (mix_Policy == "explicit") then
58 >          LJ_Mixing_Policy = Explicit_Mixing_Rule
59 >       else
60 >          write(*,*) 'Unknown Mixing Policy!'
61 >          status = -1
62 >          return
63 >       endif
64 >    endif
65 >        
66      status = 0
67 <    call createMixingList(ListHead,myStatus)
67 >    call createMixingList(myStatus)
68      if (myStatus /= 0) then
69         status = -1
70         return
71      end if
72 +    
73 +    LJ_FF_initialized = .true.
74  
75 <  end subroutine init_ljFF
76 <
77 <  subroutine createMixingList(ListHead,status)
78 <    type (atype), pointer :: ListHead
63 <    integer :: listSize
75 >  end subroutine init_LJ_FF
76 >  
77 >  subroutine createMixingList(status)
78 >    integer :: nAtypes
79      integer :: status
80      integer :: i
81      integer :: j
82 <
83 <    integer :: i = 0
84 <    integer :: j = 0
70 <    type (atype), pointer :: tmpPtr_i => null()
71 <    type (atype), pointer :: tmpPtr_j => null()
82 >    real ( kind = dp ) :: mySigma_i,mySigma_j
83 >    real ( kind = dp ) :: myEpsilon_i,myEpsilon_j
84 >    real ( kind = dp ) :: rcut, rcut6
85      status = 0
86 <
87 <    listSize = getListLen(ListHead)
88 <    if (listSize == 0) then
86 >    
87 >    nAtypes = getSize(atypes)
88 >    if (nAtypes == 0) then
89         status = -1
90         return
91      end if
92 <  
80 <
92 >        
93      if (.not. associated(ljMixed)) then
94 <       allocate(ljMixed(listSize,listSize))
94 >       allocate(ljMixed(nAtypes, nAtypes))
95      else
96         status = -1
97         return
98      end if
99 +    call getRcut(rcut, rc6=rcut6)
100 +    do i = 1, nAtypes
101  
102 <    
103 <
104 <    tmpPtr_i => ListHead
105 <    tmpPtr_j => tmpPtr_i%next
106 <    do while (associated(tmpPtr_i))
93 <       i = i + 1
94 < ! do self mixing rule
95 <       ljMixed(i,i)%sigma   = tmpPtr_i%sigma
96 <                                                                                                  
97 <       ljMixed(i,i)%sigma2  = (ljMixed(i,i)%sigma) ** 2
98 <                                                                                                  
102 >       call getElementProperty(atypes,i,"lj_epsilon",myEpsilon_i)
103 >       call getElementProperty(atypes,i,"lj_sigma",  mySigma_i)
104 >       ! do self mixing rule
105 >       ljMixed(i,i)%sigma   = mySigma_i
106 >      
107         ljMixed(i,i)%sigma6  = (ljMixed(i,i)%sigma) ** 6                                                                                                  
108 <       ljMixed(i,i)%epsilon = tmpPtr_i%epsilon
108 >       ljMixed(i,i)%tp6     = ljMixed(i,i)%sigma6/rcut6
109  
110 <       j = i + 1
111 <       do while (associated(tmpPtr_j))
110 >       ljMixed(i,i)%tp12    = (ljMixed(i,i)%tp6) ** 2
111 >
112 >       ljMixed(i,i)%epsilon = myEpsilon_i
113 >
114 >       ljMixed(i,i)%delta = -4.0_DP * ljMixed(i,i)%epsilon * &
115 >            (ljMixed(i,i)%tp12 - ljMixed(i,i)%tp6)
116 >      
117 >       do j = i + 1, nAtypes
118 >          call getElementProperty(atypes,j,"lj_epsilon",myEpsilon_j)
119 >          call getElementProperty(atypes,j,"lj_sigma",  mySigma_j)
120            
121            ljMixed(i,j)%sigma  =  &
122 <               calcLJMix("sigma",tmpPtr_i%sigma, &
123 <               tmpPtr_j%sigma)
122 >               calcLJMix("sigma",mySigma_i, &
123 >               mySigma_j)
124            
109          ljMixed(i,j)%sigma2  = &
110               (ljMixed(i,j)%sigma)**2
111          
125            ljMixed(i,j)%sigma6 = &
126                 (ljMixed(i,j)%sigma)**6
127            
128 +          
129 +          ljMixed(i,j)%tp6     = ljMixed(i,j)%sigma6/rcut6
130 +          
131 +          ljMixed(i,j)%tp12    = (ljMixed(i,j)%tp6) ** 2
132 +          
133 +          
134            ljMixed(i,j)%epsilon = &
135 <               calcLJMix("epsilon",tmpPtr_i%epsilon, &
136 <               tmpPtr_j%epsilon)
135 >               calcLJMix("epsilon",myEpsilon_i, &
136 >               myEpsilon_j)
137 >          
138 >          ljMixed(i,j)%delta = -4.0_DP * ljMixed(i,j)%epsilon * &
139 >               (ljMixed(i,j)%tp12 - ljMixed(i,j)%tp6)
140 >          
141 >          
142            ljMixed(j,i)%sigma   = ljMixed(i,j)%sigma
119          ljMixed(j,i)%sigma2  = ljMixed(i,j)%sigma2
143            ljMixed(j,i)%sigma6  = ljMixed(i,j)%sigma6
144 +          ljMixed(j,i)%tp6     = ljMixed(i,j)%tp6
145 +          ljMixed(j,i)%tp12    = ljMixed(i,j)%tp12
146            ljMixed(j,i)%epsilon = ljMixed(i,j)%epsilon
147 <
148 <
124 <          tmpPtr_j => tmpPtr_j%next
125 <          j = j + 1
147 >          ljMixed(j,i)%delta   = ljMixed(i,j)%delta
148 >          
149         end do
127 ! advance pointers
128       tmpPtr_i => tmpPtr_i%next
129       if (associated(tmpPtr_i)) then
130          tmpPtr_j => tmpPtr_i%next
131       endif
132      
150      end do
151 <
151 >    
152    end subroutine createMixingList
153 +  
154  
155 +  subroutine do_lj_pair(atom1, atom2, d, rij, r2, pot, f, do_pot, do_stress)
156  
157 <
158 <
157 >    integer, intent(in) ::  atom1, atom2
158 >    real( kind = dp ), intent(in) :: rij, r2
159 >    real( kind = dp ) :: pot
160 >    real( kind = dp ), dimension(3, getNlocal()) :: f    
161 >    real( kind = dp ), intent(in), dimension(3) :: d
162 >    logical, intent(in) :: do_pot, do_stress
163  
164 < !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
165 < !! derivatives.
166 <  subroutine getLjForce(r,pot,dudr,atype1,atype2,d2,status)
167 < ! arguments
145 < !! Length of vector between particles
146 <    real( kind = dp ), intent(in)  :: r
147 < !! Potential Energy
148 <    real( kind = dp ), intent(out) :: pot
149 < !! Derivatve wrt postion
150 <    real( kind = dp ), intent(out) :: dudr
151 < !! Second Derivative, optional, used mainly for normal mode calculations.
152 <    real( kind = dp ), intent(out), optional :: d2
153 <    
154 <    type (atype), pointer :: atype1
155 <    type (atype), pointer :: atype2
156 <
157 <    integer, intent(out), optional :: status
158 <
159 < ! local Variables
160 <    real( kind = dp ) :: sigma
161 <    real( kind = dp ) :: sigma2
164 >    ! local Variables
165 >    real( kind = dp ) :: drdx, drdy, drdz
166 >    real( kind = dp ) :: fx, fy, fz
167 >    real( kind = dp ) :: pot_temp, dudr
168      real( kind = dp ) :: sigma6
169      real( kind = dp ) :: epsilon
164
170      real( kind = dp ) :: rcut
166    real( kind = dp ) :: rcut2
167    real( kind = dp ) :: rcut6
168    real( kind = dp ) :: r2
171      real( kind = dp ) :: r6
170
172      real( kind = dp ) :: t6
173      real( kind = dp ) :: t12
173    real( kind = dp ) :: tp6
174    real( kind = dp ) :: tp12
174      real( kind = dp ) :: delta
175  
176 <    logical :: doSec = .false.
177 <
178 <    integer :: errorStat
179 <
180 < !! Optional Argument Checking
181 < ! Check to see if we need to do second derivatives
176 >    ! Look up the correct parameters in the mixing matrix
177 > #ifdef IS_MPI
178 >    sigma6   = ljMixed(atid_Row(atom1),atid_Col(atom2))%sigma6
179 >    epsilon  = ljMixed(atid_Row(atom1),atid_Col(atom2))%epsilon
180 >    delta    = ljMixed(atid_Row(atom1),atid_Col(atom2))%delta
181 > #else
182 >    sigma6   = ljMixed(atid(atom1),atid(atom2))%sigma6
183 >    epsilon  = ljMixed(atid(atom1),atid(atom2))%epsilon
184 >    delta    = ljMixed(atid(atom1),atid(atom2))%delta
185 > #endif
186 >  
187 >    call getRcut(rcut)
188      
184    if (present(d2))     doSec = .true.
185    if (present(status)) status = 0
186
187 ! Look up the correct parameters in the mixing matrix
188    sigma    = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma
189    sigma2   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma2
190    sigma6   = ljMixed(atype1%atype_ident,atype2%atype_ident)%sigma6
191    epsilon  = ljMixed(atype1%atype_ident,atype2%atype_ident)%epsilon
192
193
194    
195
196    call getRcut(rcut,rcut2=rcut2,rcut6=rcut6,status=errorStat)
197    
198    r2 = r * r
189      r6 = r2 * r2 * r2
190  
191      t6  = sigma6/ r6
192      t12 = t6 * t6
193 +      
194 +    if (rij.le.rcut) then
195 +      
196 +       pot_temp = 4.0E0_DP * epsilon * (t12 - t6) + delta      
197  
198 +       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
199  
200 <                                                                              
201 <    tp6 = sigma6 / rcut6
202 <    tp12 = tp6*tp6
203 <                                                                              
204 <    delta = -4.0E0_DP*epsilon * (tp12 - tp6)
205 <                                                                              
206 <    if (r.le.rcut) then
207 <       pot = 4.0E0_DP * epsilon * (t12 - t6) + delta
208 <       dudr = 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / r
209 <       if(doSec)  d2 = 24.0E0_DP * epsilon * (26.0E0_DP*t12 - 7.0E0_DP*t6)/r/r
210 <    else
211 <       pot = 0.0E0_DP
212 <       dudr = 0.0E0_DP
218 <       if(doSec) d2 = 0.0E0_DP
219 <    endif
220 <                                                                              
221 <  return
200 >       drdx = -d(1) / rij
201 >       drdy = -d(2) / rij
202 >       drdz = -d(3) / rij
203 >      
204 >       fx = dudr * drdx
205 >       fy = dudr * drdy
206 >       fz = dudr * drdz
207 >    
208 > #ifdef IS_MPI
209 >       if (do_pot) then
210 >          pot_Row(atom1) = pot_Row(atom1) + pot_temp*0.5
211 >          pot_Col(atom2) = pot_Col(atom2) + pot_temp*0.5
212 >       endif
213  
214 +       f_Row(1,atom1) = f_Row(1,atom1) + fx
215 +       f_Row(2,atom1) = f_Row(2,atom1) + fy
216 +       f_Row(3,atom1) = f_Row(3,atom1) + fz
217  
218 +       f_Col(1,atom2) = f_Col(1,atom2) - fx
219 +       f_Col(2,atom2) = f_Col(2,atom2) - fy
220 +       f_Col(3,atom2) = f_Col(3,atom2) - fz      
221  
222 <  end subroutine getLJForce
222 > #else
223 >       if (do_pot) pot = pot + pot_temp
224  
225 +       f(1,atom1) = f(1,atom1) + fx
226 +       f(2,atom1) = f(2,atom1) + fy
227 +       f(3,atom1) = f(3,atom1) + fz
228  
229 < !! Calculates the mixing for sigma or epslon based on character string for initialzition of mixing array.
229 >       f(1,atom2) = f(1,atom2) - fx
230 >       f(2,atom2) = f(2,atom2) - fy
231 >       f(3,atom2) = f(3,atom2) - fz
232 > #endif
233 >      
234 >       if (do_stress) then
235 >          tau_Temp(1) = tau_Temp(1) + fx * d(1)
236 >          tau_Temp(2) = tau_Temp(2) + fx * d(2)
237 >          tau_Temp(3) = tau_Temp(3) + fx * d(3)
238 >          tau_Temp(4) = tau_Temp(4) + fy * d(1)
239 >          tau_Temp(5) = tau_Temp(5) + fy * d(2)
240 >          tau_Temp(6) = tau_Temp(6) + fy * d(3)
241 >          tau_Temp(7) = tau_Temp(7) + fz * d(1)
242 >          tau_Temp(8) = tau_Temp(8) + fz * d(2)
243 >          tau_Temp(9) = tau_Temp(9) + fz * d(3)
244 >          virial_Temp = virial_Temp + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
245 >       endif
246 >      
247 >    endif
248 >    return          
249 >  end subroutine do_lj_pair
250 >
251 >
252 >  !! Calculates the mixing for sigma or epslon
253 >
254    function calcLJMix(thisParam,param1,param2,status) result(myMixParam)
255      character(len=*) :: thisParam
256      real(kind = dp)  :: param1
257      real(kind = dp)  :: param2
258      real(kind = dp ) :: myMixParam
234    character(len = getStringLen()) :: thisMixingRule
235    integer, optional :: status
259  
260 < !! get the mixing rules from the simulation
238 <    thisMixingRule = returnMixingRules()
239 <    myMixParam = 0.0_dp
260 >    integer, optional :: status  
261  
262 +    myMixParam = 0.0_dp
263 +    
264      if (present(status)) status = 0
265 <    select case (thisMixingRule)
266 <       case ("standard")
267 <          select case (thisParam)
268 <          case ("sigma")
269 <             myMixParam = 0.5_dp * (param1 + param2)
270 <          case ("epsilon")
271 <             myMixParam = sqrt(param1 * param2)
272 <          case default
273 <             status = -1
274 <          end select
252 <    case("LJglass")
265 >    select case (LJ_Mixing_Policy)
266 >    case (LB_Mixing_Rule)
267 >       select case (thisParam)
268 >       case ("sigma")
269 >          myMixParam = 0.5_dp * (param1 + param2)
270 >       case ("epsilon")
271 >          myMixParam = sqrt(param1 * param2)
272 >       case default
273 >          status = -1
274 >       end select
275      case default
276         status = -1
277      end select
278    end function calcLJMix
279 <
280 <
259 <
260 < end module lj_ff
279 >  
280 > end module lj

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines