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 1948 by gezelter, Fri Jan 14 20:31:16 2005 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 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.7 2005-01-14 20:31:16 gezelter Exp $, $Date: 2005-01-14 20:31:16 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $
46 > !! @version $Id: LJ.F90,v 1.11 2005-04-15 22:03:47 gezelter Exp $, $Date: 2005-04-15 22:03:47 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $
47  
48  
49   module lj
# Line 59 | Line 59 | module lj
59  
60    implicit none
61    PRIVATE
62 <  
62 >
63    integer, parameter :: DP = selected_real_kind(15)
64 <  
64 >
65    type, private :: LjType
66       integer :: c_ident
67       integer :: atid
68       real(kind=dp) :: sigma
69       real(kind=dp) :: epsilon
70 +     logical :: soft_pot
71    end type LjType
72 <  
72 >
73    type(LjType), dimension(:), allocatable :: ParameterMap
74 <  
74 >
75    logical, save :: haveMixingMap = .false.
76 <  
76 >
77    type :: MixParameters
78       real(kind=DP) :: sigma
79       real(kind=DP) :: epsilon
# Line 80 | Line 81 | module lj
81       real(kind=dp)  :: tp6
82       real(kind=dp)  :: tp12
83       real(kind=dp)  :: delta
84 +     logical :: soft_pot    
85    end type MixParameters
86 <  
86 >
87    type(MixParameters), dimension(:,:), allocatable :: MixingMap
88 <  
88 >
89    real(kind=DP), save :: LJ_rcut
90    logical, save :: have_rcut = .false.
91    logical, save :: LJ_do_shift = .false.
92    logical, save :: useGeometricDistanceMixing = .false.
93 <  
93 >
94    !! Public methods and data
95 <  
95 >
96    public :: setCutoffLJ
97    public :: useGeometricMixing
98    public :: do_lj_pair
99    public :: newLJtype  
100    public :: getSigma
101    public :: getEpsilon
102 <  
102 >  public :: destroyLJTypes
103 >
104   contains
105  
106 <  subroutine newLJtype(c_ident, sigma, epsilon, status)
106 >  subroutine newLJtype(c_ident, sigma, epsilon, soft_pot, status)
107      integer,intent(in) :: c_ident
108      real(kind=dp),intent(in) :: sigma
109      real(kind=dp),intent(in) :: epsilon
110 +    integer, intent(in) :: soft_pot
111      integer,intent(out) :: status
112      integer :: nATypes, myATID
113      integer, pointer :: MatchList(:) => null()
114  
115      status = 0
116 <    
116 >
117      myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
118  
119      if (.not.allocated(ParameterMap)) then
120 <      
120 >
121         !call getMatchingElementList(atypes, "is_LennardJones", .true., &
122         !     nLJTypes, MatchList)
123         nAtypes = getSize(atypes)
# Line 121 | Line 125 | contains
125            status = -1
126            return
127         end if
128 <      
128 >
129         if (.not. allocated(ParameterMap)) then
130            allocate(ParameterMap(nAtypes))
131         endif
132 <      
132 >
133      end if
134  
135      if (myATID .gt. size(ParameterMap)) then
136         status = -1
137         return
138      endif
139 <    
139 >
140      ! set the values for ParameterMap for this atom type:
141  
142      ParameterMap(myATID)%c_ident = c_ident
143      ParameterMap(myATID)%atid = myATID
144      ParameterMap(myATID)%epsilon = epsilon
145      ParameterMap(myATID)%sigma = sigma
146 <    
146 >    if (soft_pot == 1) then
147 >       ParameterMap(myATID)%soft_pot = .true.
148 >    else
149 >       ParameterMap(myATID)%soft_pot = .false.
150 >    endif
151    end subroutine newLJtype
152  
153    function getSigma(atid) result (s)
154      integer, intent(in) :: atid
155      integer :: localError
156      real(kind=dp) :: s
157 <    
157 >
158      if (.not.allocated(ParameterMap)) then
159         call handleError("LJ", "no ParameterMap was present before first call of getSigma!")
160         return
161      end if
162 <    
162 >
163      s = ParameterMap(atid)%sigma
164    end function getSigma
165  
# Line 159 | Line 167 | contains
167      integer, intent(in) :: atid
168      integer :: localError
169      real(kind=dp) :: e
170 <    
170 >
171      if (.not.allocated(ParameterMap)) then
172         call handleError("LJ", "no ParameterMap was present before first call of getEpsilon!")
173         return
174      end if
175 <    
175 >
176      e = ParameterMap(atid)%epsilon
177    end function getEpsilon
178  
# Line 183 | Line 191 | contains
191      LJ_do_shift = do_shift
192      call set_switch(LJ_SWITCH, rcut, rcut)
193      have_rcut = .true.
194 <    
194 >
195      return
196    end subroutine setCutoffLJ
197  
# Line 192 | Line 200 | contains
200      haveMixingMap = .false.
201      return
202    end subroutine useGeometricMixing
203 <  
203 >
204    subroutine createMixingMap(status)
205      integer :: nATIDs
206      integer :: status
# Line 201 | Line 209 | contains
209      real ( kind = dp ) :: Sigma_i, Sigma_j
210      real ( kind = dp ) :: Epsilon_i, Epsilon_j
211      real ( kind = dp ) :: rcut6
212 +    logical :: i_is_LJ, j_is_LJ
213  
214      status = 0
215 <    
215 >
216 >    if (.not. allocated(ParameterMap)) then
217 >       call handleError("LJ", "no ParameterMap was present before call of createMixingMap!")
218 >       status = -1
219 >       return
220 >    endif
221 >
222      nATIDs = size(ParameterMap)
223 <    
223 >
224      if (nATIDs == 0) then
225         status = -1
226         return
227      end if
228  
229 +    if (.not. allocated(MixingMap)) then
230 +       allocate(MixingMap(nATIDs, nATIDs))
231 +    endif
232 +
233      if (.not.have_rcut) then
234         status = -1
235         return
236      endif
237 <    
219 <    if (.not. allocated(MixingMap)) then
220 <       allocate(MixingMap(nATIDs, nATIDs))
221 <    endif
222 <    
237 >
238      rcut6 = LJ_rcut**6
239 <    
239 >
240      ! This loops through all atypes, even those that don't support LJ forces.
241      do i = 1, nATIDs
242 <      
243 <       Epsilon_i = ParameterMap(i)%epsilon
244 <       Sigma_i = ParameterMap(i)%sigma
245 <      
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 <      
255 <       do j = i + 1, nATIDs
256 <          
257 <          Epsilon_j = ParameterMap(j)%epsilon
258 <          Sigma_j = ParameterMap(j)%sigma
259 <          
260 <          ! only the distance parameter uses different mixing policies
261 <          if (useGeometricDistanceMixing) then
262 <             ! only for OPLS as far as we can tell
263 <             MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
264 <          else
265 <             ! everyone else
266 <             MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
267 <          endif
268 <          
269 <          ! energy parameter is always geometric mean:
270 <          MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
271 <                    
272 <          MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
273 <          MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
274 <          MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
275 <          
276 <          MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
277 <               (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
278 <          
279 <          MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
280 <          MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
281 <          MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
282 <          MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
283 <          MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
284 <          MixingMap(j,i)%delta   = MixingMap(i,j)%delta
285 <          
286 <       end do
242 >       call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
243 >       if (i_is_LJ) then
244 >          Epsilon_i = ParameterMap(i)%epsilon
245 >          Sigma_i = ParameterMap(i)%sigma
246 >
247 >          ! do self mixing rule
248 >          MixingMap(i,i)%sigma   = Sigma_i          
249 >          MixingMap(i,i)%sigma6  = Sigma_i ** 6          
250 >          MixingMap(i,i)%tp6     = (MixingMap(i,i)%sigma6)/rcut6          
251 >          MixingMap(i,i)%tp12    = (MixingMap(i,i)%tp6) ** 2
252 >          MixingMap(i,i)%epsilon = Epsilon_i          
253 >          MixingMap(i,i)%delta   = -4.0_DP * MixingMap(i,i)%epsilon * &
254 >               (MixingMap(i,i)%tp12 - MixingMap(i,i)%tp6)
255 >          MixingMap(i,i)%soft_pot = ParameterMap(i)%soft_pot
256 >
257 >          do j = i + 1, nATIDs
258 >             call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
259 >
260 >             if (j_is_LJ) then
261 >                Epsilon_j = ParameterMap(j)%epsilon
262 >                Sigma_j = ParameterMap(j)%sigma
263 >
264 >                ! only the distance parameter uses different mixing policies
265 >                if (useGeometricDistanceMixing) then
266 >                   ! only for OPLS as far as we can tell
267 >                   MixingMap(i,j)%sigma = dsqrt(Sigma_i * Sigma_j)
268 >                else
269 >                   ! everyone else
270 >                   MixingMap(i,j)%sigma = 0.5_dp * (Sigma_i + Sigma_j)
271 >                endif
272 >
273 >                ! energy parameter is always geometric mean:
274 >                MixingMap(i,j)%epsilon = dsqrt(Epsilon_i * Epsilon_j)
275 >
276 >                MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
277 >                MixingMap(i,j)%tp6    = MixingMap(i,j)%sigma6/rcut6
278 >                MixingMap(i,j)%tp12    = (MixingMap(i,j)%tp6) ** 2
279 >
280 >                MixingMap(i,j)%delta = -4.0_DP * MixingMap(i,j)%epsilon * &
281 >                     (MixingMap(i,j)%tp12 - MixingMap(i,j)%tp6)
282 >
283 >                MixingMap(i,j)%soft_pot = ParameterMap(i)%soft_pot .or. ParameterMap(j)%soft_pot
284 >
285 >
286 >                MixingMap(j,i)%sigma   = MixingMap(i,j)%sigma
287 >                MixingMap(j,i)%sigma6  = MixingMap(i,j)%sigma6
288 >                MixingMap(j,i)%tp6     = MixingMap(i,j)%tp6
289 >                MixingMap(j,i)%tp12    = MixingMap(i,j)%tp12
290 >                MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon
291 >                MixingMap(j,i)%delta   = MixingMap(i,j)%delta
292 >                MixingMap(j,i)%soft_pot   = MixingMap(i,j)%soft_pot
293 >             endif
294 >          end do
295 >       endif
296      end do
297 <    
297 >
298      haveMixingMap = .true.
299  
300    end subroutine createMixingMap
301 <        
301 >
302    subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
303         pot, f, do_pot)
304  
# Line 296 | Line 320 | contains
320      real( kind = dp ) :: t6
321      real( kind = dp ) :: t12
322      real( kind = dp ) :: delta
323 +    logical :: soft_pot
324      integer :: id1, id2, localError
325  
326      if (.not.haveMixingMap) then
# Line 312 | Line 337 | contains
337      sigma6   = MixingMap(atid_Row(atom1),atid_Col(atom2))%sigma6
338      epsilon  = MixingMap(atid_Row(atom1),atid_Col(atom2))%epsilon
339      delta    = MixingMap(atid_Row(atom1),atid_Col(atom2))%delta
340 +    soft_pot =  MixingMap(atid_Row(atom1),atid_Col(atom2))%soft_pot
341   #else
342      sigma6   = MixingMap(atid(atom1),atid(atom2))%sigma6
343      epsilon  = MixingMap(atid(atom1),atid(atom2))%epsilon
344      delta    = MixingMap(atid(atom1),atid(atom2))%delta
345 +    soft_pot    = MixingMap(atid(atom1),atid(atom2))%soft_pot
346   #endif
347  
348      r6 = r2 * r2 * r2
349 <    
349 >
350      t6  = sigma6/ r6
351      t12 = t6 * t6    
352 <  
353 <    pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
354 <    if (LJ_do_shift) then
355 <       pot_temp = pot_temp + delta
352 >
353 >    if (soft_pot) then
354 >       pot_temp = 4.0E0_DP * epsilon * t12
355 >       if (LJ_do_shift) then
356 >          pot_temp = pot_temp + delta
357 >       endif
358 >
359 >       vpair = vpair + pot_temp
360 >
361 >       dudr = sw * 24.0E0_DP * epsilon * (-2.0E0_DP)*t12 / rij
362 >
363 >    else
364 >       pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
365 >       if (LJ_do_shift) then
366 >          pot_temp = pot_temp + delta
367 >       endif
368 >
369 >       vpair = vpair + pot_temp
370 >
371 >       dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
372      endif
373  
331 –    vpair = vpair + pot_temp
332 –      
333 –    dudr = sw * 24.0E0_DP * epsilon * (t6 - 2.0E0_DP*t12) / rij
334 –      
374      drdx = d(1) / rij
375      drdy = d(2) / rij
376      drdz = d(3) / rij
377 <      
377 >
378      fx = dudr * drdx
379      fy = dudr * drdy
380      fz = dudr * drdz
381 <    
382 <      
381 >
382 >
383   #ifdef IS_MPI
384      if (do_pot) then
385         pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
386         pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
387      endif
388 <    
388 >
389      f_Row(1,atom1) = f_Row(1,atom1) + fx
390      f_Row(2,atom1) = f_Row(2,atom1) + fy
391      f_Row(3,atom1) = f_Row(3,atom1) + fz
392 <    
392 >
393      f_Col(1,atom2) = f_Col(1,atom2) - fx
394      f_Col(2,atom2) = f_Col(2,atom2) - fy
395      f_Col(3,atom2) = f_Col(3,atom2) - fz      
396 <    
396 >
397   #else
398      if (do_pot) pot = pot + sw*pot_temp
399  
400      f(1,atom1) = f(1,atom1) + fx
401      f(2,atom1) = f(2,atom1) + fy
402      f(3,atom1) = f(3,atom1) + fz
403 <    
403 >
404      f(1,atom2) = f(1,atom2) - fx
405      f(2,atom2) = f(2,atom2) - fy
406      f(3,atom2) = f(3,atom2) - fz
407   #endif
408 <        
408 >
409   #ifdef IS_MPI
410      id1 = AtomRowToGlobal(atom1)
411      id2 = AtomColToGlobal(atom2)
# Line 376 | Line 415 | contains
415   #endif
416  
417      if (molMembershipList(id1) .ne. molMembershipList(id2)) then
418 <      
418 >
419         fpair(1) = fpair(1) + fx
420         fpair(2) = fpair(2) + fy
421         fpair(3) = fpair(3) + fz
# Line 384 | Line 423 | contains
423      endif
424  
425      return    
426 <    
426 >
427    end subroutine do_lj_pair
428 <  
429 <  
430 <  !! Calculates the mixing for sigma or epslon
431 <    
428 >
429 >  subroutine destroyLJTypes()
430 >    if(allocated(ParameterMap)) deallocate(ParameterMap)
431 >    if(allocated(MixingMap)) deallocate(MixingMap)
432 >    haveMixingMap = .false.
433 >  end subroutine destroyLJTypes
434 >
435 >
436   end module lj

Diff Legend

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