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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines