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 2085 by tim, Fri Feb 25 21:22:00 2005 UTC vs.
Revision 2086 by gezelter, Tue Mar 8 21:06:12 2005 UTC

# Line 43 | Line 43
43   !! Calculates Long Range forces Lennard-Jones interactions.
44   !! @author Charles F. Vardeman II
45   !! @author Matthew Meineke
46 < !! @version $Id: LJ.F90,v 1.8 2005-02-25 21:22:00 tim Exp $, $Date: 2005-02-25 21:22:00 $, $Name: not supported by cvs2svn $, $Revision: 1.8 $
46 > !! @version $Id: LJ.F90,v 1.9 2005-03-08 21:06:12 gezelter Exp $, $Date: 2005-03-08 21:06:12 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $
47  
48  
49   module lj
# Line 124 | Line 124 | contains
124            status = -1
125            return
126         end if
127 <      
127 >        
128         if (.not. allocated(ParameterMap)) then
129            allocate(ParameterMap(nAtypes))
130         endif
131 <      
131 >      
132      end if
133  
134      if (myATID .gt. size(ParameterMap)) then
# Line 208 | Line 208 | contains
208      real ( kind = dp ) :: Sigma_i, Sigma_j
209      real ( kind = dp ) :: Epsilon_i, Epsilon_j
210      real ( kind = dp ) :: rcut6
211 +    logical :: i_is_LJ, j_is_LJ
212  
213      status = 0
214 +
215 +    if (.not. allocated(ParameterMap)) then
216 +       call handleError("LJ", "no ParameterMap was present before call of createMixingMap!")
217 +       status = -1
218 +       return
219 +    endif
220      
221      nATIDs = size(ParameterMap)
222      
# Line 218 | Line 225 | contains
225         return
226      end if
227  
228 +    if (.not. allocated(MixingMap)) then
229 +       allocate(MixingMap(nATIDs, nATIDs))
230 +    endif
231 +
232      if (.not.have_rcut) then
233         status = -1
234         return
235      endif
236 <    
226 <    if (.not. allocated(MixingMap)) then
227 <       allocate(MixingMap(nATIDs, nATIDs))
228 <    endif
229 <    
236 >        
237      rcut6 = LJ_rcut**6
238      
239      ! This loops through all atypes, even those that don't support LJ forces.
240      do i = 1, nATIDs
241 <      
242 <       Epsilon_i = ParameterMap(i)%epsilon
243 <       Sigma_i = ParameterMap(i)%sigma
244 <      
238 <       ! do self mixing rule
239 <       MixingMap(i,i)%sigma   = Sigma_i          
240 <       MixingMap(i,i)%sigma6  = Sigma_i ** 6          
241 <       MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
242 <       MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
243 <       MixingMap(i,i)%epsilon = Epsilon_i          
244 <       MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
245 <            (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
246 <       MixingMap(i,i)%soft_pot = ParameterMap(i)%soft_pot
247 <      
248 <       do j = i + 1, nATIDs
241 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
242 >       if (i_is_LJ) then
243 >          Epsilon_i = ParameterMap(i)%epsilon
244 >          Sigma_i = ParameterMap(i)%sigma
245            
246 <          Epsilon_j = ParameterMap(j)%epsilon
247 <          Sigma_j = ParameterMap(j)%sigma
246 >          ! do self mixing rule
247 >          MixingMap(i,i)%sigma   = Sigma_i          
248 >          MixingMap(i,i)%sigma6  = Sigma_i ** 6          
249 >          MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
250 >          MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
251 >          MixingMap(i,i)%epsilon = Epsilon_i          
252 >          MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
253 >               (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
254 >          MixingMap(i,i)%soft_pot = ParameterMap(i)%soft_pot
255            
256 <          ! only the distance parameter uses different mixing policies
257 <          if (useGeometricDistanceMixing) then
255 <             ! only for OPLS as far as we can tell
256 <             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
257 <          else
258 <             ! everyone else
259 <             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
260 <          endif
256 >          do j = i + 1, nATIDs
257 >             call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
258            
259 <          ! energy parameter is always geometric mean:
260 <          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
261 <                    
262 <          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
263 <          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
264 <          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
265 <          
266 <          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
267 <               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
268 <
269 <          MixingMap(i,j)%soft_pot = ParameterMap(i)%soft_pot .or. ParameterMap(j)%soft_pot
270 <
271 <          
272 <          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
273 <          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
274 <          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
275 <          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
276 <          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
277 <          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
278 <          MixingMap(j,i)%soft_pot   = MixingMap(i,j)%soft_pot
279 <          
280 <       end do
259 >             if (j_is_LJ) then
260 >                Epsilon_j = ParameterMap(j)%epsilon
261 >                Sigma_j = ParameterMap(j)%sigma
262 >                
263 >                ! only the distance parameter uses different mixing policies
264 >                if (useGeometricDistanceMixing) then
265 >                   ! only for OPLS as far as we can tell
266 >                   MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
267 >                else
268 >                   ! everyone else
269 >                   MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
270 >                endif
271 >                
272 >                ! energy parameter is always geometric mean:
273 >                MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
274 >                
275 >                MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
276 >                MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
277 >                MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
278 >                
279 >                MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
280 >                     (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
281 >                
282 >                MixingMap(i,j)%soft_pot = ParameterMap(i)%soft_pot .or. ParameterMap(j)%soft_pot
283 >                
284 >                
285 >                MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
286 >                MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
287 >                MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
288 >                MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
289 >                MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
290 >                MixingMap(j,i)%delta   = MixingMap(i,j)%delta
291 >                MixingMap(j,i)%soft_pot   = MixingMap(i,j)%soft_pot
292 >             endif
293 >          end do
294 >       endif
295      end do
296      
297      haveMixingMap = .true.

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines