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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines