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 2319 by chuckv, Wed Sep 21 23:45:48 2005 UTC vs.
Revision 2461 by gezelter, Mon Nov 21 22:59:02 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.16 2005-09-21 23:45:48 chuckv Exp $, $Date: 2005-09-21 23:45:48 $, $Name: not supported by cvs2svn $, $Revision: 1.16 $
46 > !! @version $Id: LJ.F90,v 1.19 2005-11-21 22:59:01 gezelter Exp $, $Date: 2005-11-21 22:59:01 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
47  
48  
49   module lj
# Line 58 | Line 58 | module lj
58  
59    implicit none
60    PRIVATE
61 + #define __FORTRAN90
62 + #include "UseTheForce/DarkSide/fInteractionMap.h"
63  
64    integer, parameter :: DP = selected_real_kind(15)
65  
# Line 100 | Line 102 | module lj
102  
103    public :: newLJtype
104    public :: setLJDefaultCutoff
103 –  public :: setLJUniformCutoff
104 –  public :: setLJCutoffByTypes
105    public :: getSigma
106    public :: getEpsilon
107    public :: useGeometricMixing
# Line 162 | Line 162 | contains
162         call createMixingMap()
163      end if
164    end subroutine setLJDefaultCutoff
165 –
166 –  subroutine setLJUniformCutoff(thisRcut, shiftedPot)
167 –    real(kind=dp), intent(in) :: thisRcut
168 –    logical, intent(in) :: shiftedPot
169 –    integer :: nLJtypes, i, j
170 –
171 –    if (LJMap%currentLJtype == 0) then
172 –       call handleError("LJ", "No members in LJMap")
173 –       return
174 –    end if
175 –
176 –    nLJtypes = LJMap%nLJtypes
177 –    if (.not. allocated(MixingMap)) then
178 –       allocate(MixingMap(nLJtypes, nLJtypes))
179 –    endif
180 –
181 –    do i = 1, nLJtypes
182 –       do j = 1, nLJtypes
183 –          MixingMap(i,j)%rCut = thisRcut
184 –          MixingMap(i,j)%shiftedPot = shiftedPot
185 –          MixingMap(i,j)%rCutWasSet = .true.
186 –       enddo
187 –    enddo
188 –    call createMixingMap()
189 –  end subroutine setLJUniformCutoff
190 –
191 –  subroutine setLJCutoffByTypes(atid1, atid2, thisRcut, shiftedPot)
192 –    integer, intent(in) :: atid1, atid2
193 –    real(kind=dp), intent(in) :: thisRcut
194 –    logical, intent(in) :: shiftedPot
195 –    integer :: nLJtypes, ljt1, ljt2
196 –
197 –    if (LJMap%currentLJtype == 0) then
198 –       call handleError("LJ", "No members in LJMap")
199 –       return
200 –    end if
201 –
202 –    nLJtypes = LJMap%nLJtypes
203 –    if (.not. allocated(MixingMap)) then
204 –       allocate(MixingMap(nLJtypes, nLJtypes))
205 –    endif
206 –
207 –    ljt1 = LJMap%atidToLJtype(atid1)
208 –    ljt2 = LJMap%atidToLJtype(atid2)
209 –
210 –    MixingMap(ljt1,ljt2)%rCut = thisRcut
211 –    MixingMap(ljt1,ljt2)%shiftedPot = shiftedPot
212 –    MixingMap(ljt1,ljt2)%rCutWasSet = .true.
213 –    MixingMap(ljt2,ljt1)%rCut = thisRcut
214 –    MixingMap(ljt2,ljt1)%shiftedPot = shiftedPot
215 –    MixingMap(ljt2,ljt1)%rCutWasSet = .true.
165  
217 –    call createMixingMap()
218 –  end subroutine setLJCutoffByTypes
219 –
166    function getSigma(atid) result (s)
167      integer, intent(in) :: atid
168      integer :: ljt1
# Line 295 | Line 241 | contains
241  
242            MixingMap(i,j)%sigma6 = (MixingMap(i,j)%sigma)**6
243  
244 <          if (MixingMap(i,j)%rCutWasSet) then
245 <             rcut6 = (MixingMap(i,j)%rcut)**6
244 >          if (haveDefaultCutoff) then
245 >             rcut6 = defaultCutoff**6
246 >             tp6    = MixingMap(i,j)%sigma6/rcut6
247 >             tp12    = tp6**2          
248 >             MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
249 >             MixingMap(i,j)%shiftedPot = defaultShift
250            else
251 <             if (haveDefaultCutoff) then
252 <                rcut6 = defaultCutoff**6
253 <                doShift = defaultShift
304 <             else
305 <                call handleError("LJ", "No specified or default cutoff value!")
306 <             endif
307 <          endif
308 <          
309 <          tp6    = MixingMap(i,j)%sigma6/rcut6
310 <          tp12    = tp6**2          
311 <          MixingMap(i,j)%delta =-4.0_DP*MixingMap(i,j)%epsilon*(tp12 - tp6)
312 <          MixingMap(i,j)%shiftedPot = doShift
251 >             MixingMap(i,j)%delta = 0.0_DP
252 >             MixingMap(i,j)%shiftedPot = defaultShift
253 >          endif          
254  
255            if (i.ne.j) then
256               MixingMap(j,i)%sigma      = MixingMap(i,j)%sigma
# Line 329 | Line 270 | contains
270      
271    end subroutine createMixingMap
272    
273 <  subroutine do_lj_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
273 >  subroutine do_lj_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
274         pot, f, do_pot)
275 <
275 >    
276      integer, intent(in) ::  atom1, atom2
277      integer :: atid1, atid2, ljt1, ljt2
278 <    real( kind = dp ), intent(in) :: rij, r2
278 >    real( kind = dp ), intent(in) :: rij, r2, rcut
279      real( kind = dp ) :: pot, sw, vpair
280      real( kind = dp ), dimension(3,nLocal) :: f    
281      real( kind = dp ), intent(in), dimension(3) :: d
# Line 347 | Line 288 | contains
288      real( kind = dp ) :: pot_temp, dudr
289      real( kind = dp ) :: sigma6
290      real( kind = dp ) :: epsilon
291 <    real( kind = dp ) :: r6
292 <    real( kind = dp ) :: t6
293 <    real( kind = dp ) :: t12
291 >    real( kind = dp ) :: r6, rc6
292 >    real( kind = dp ) :: t6, tp6
293 >    real( kind = dp ) :: t12, tp12
294      real( kind = dp ) :: delta
295      logical :: isSoftCore, shiftedPot
296      integer :: id1, id2, localError
# Line 385 | Line 326 | contains
326        
327         pot_temp = 4.0E0_DP * epsilon * t6
328         if (shiftedPot) then
329 +          rc6 = rcut**6
330 +          tp6 = sigma6 / rc6
331 +          delta =-4.0_DP*epsilon*(tp6)
332            pot_temp = pot_temp + delta
333         endif
334        
# Line 395 | Line 339 | contains
339      else
340         pot_temp = 4.0E0_DP * epsilon * (t12 - t6)
341         if (shiftedPot) then
342 +          rc6 = rcut**6
343 +          tp6 = sigma6 / rc6
344 +          tp12 = tp6*tp6
345 +          delta =-4.0_DP*epsilon*(tp12 - tp6)
346            pot_temp = pot_temp + delta
347         endif
348        
# Line 413 | Line 361 | contains
361  
362   #ifdef IS_MPI
363      if (do_pot) then
364 <       pot_Row(atom1) = pot_Row(atom1) + sw*pot_temp*0.5
365 <       pot_Col(atom2) = pot_Col(atom2) + sw*pot_temp*0.5
364 >       pot_Row(VDW_POT,atom1) = pot_Row(VDW_POT,atom1) + sw*pot_temp*0.5
365 >       pot_Col(VDW_POT,atom2) = pot_Col(VDW_POT,atom2) + sw*pot_temp*0.5
366      endif
367  
368      f_Row(1,atom1) = f_Row(1,atom1) + fx

Diff Legend

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