ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/suttonchen.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/suttonchen.F90 (file contents):
Revision 2430 by chuckv, Mon Nov 14 22:20:35 2005 UTC vs.
Revision 2541 by chuckv, Tue Jan 10 21:14:32 2006 UTC

# Line 1 | Line 1
1 < !!
1 > !!
2   !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   !!
4   !! The University of Notre Dame grants you ("Licensee") a
# Line 49 | Line 49 | module suttonchen
49    use status
50    use atype_module
51    use vector_class
52 +  use fForceOptions
53   #ifdef IS_MPI
54    use mpiSimulation
55   #endif
# Line 65 | Line 66 | module suttonchen
66    logical, save :: haveRcut = .false.
67    logical, save :: haveMixingMap = .false.
68    logical, save :: useGeometricDistanceMixing = .false.
69 +  logical, save :: cleanArrays = .true.
70 +  logical, save :: arraysAllocated = .false.
71  
72  
70 –
71 –
73    character(len = statusMsgSize) :: errMesg
74    integer :: sc_err
75  
76    character(len = 200) :: errMsg
77    character(len=*), parameter :: RoutineName =  "Sutton-Chen MODULE"
78    !! Logical that determines if eam arrays should be zeroed
78 –  logical :: cleanme = .true.
79    logical :: nmflag  = .false.
80  
81  
# Line 140 | Line 140 | module suttonchen
140    public :: do_SC_pair
141    public :: newSCtype
142    public :: calc_SC_prepair_rho
143 +  public :: calc_SC_preforce_Frho
144    public :: clean_SC
145    public :: destroySCtypes
146    public :: getSCCut
147   ! public :: setSCDefaultCutoff
148   ! public :: setSCUniformCutoff
149 <  public :: useGeometricMixing
149 >
150  
151   contains
152  
# Line 175 | Line 176 | contains
176  
177      ! check to see if this is the first time into
178      if (.not.associated(SCList%SCTypes)) then
179 <       call getMatchingElementList(atypes, "is_SuttonChen", .true., nSCtypes, MatchList)
179 >       call getMatchingElementList(atypes, "is_SC", .true., nSCtypes, MatchList)
180         SCList%nSCtypes = nSCtypes
181         allocate(SCList%SCTypes(nSCTypes))
182         nAtypes = getSize(atypes)
# Line 220 | Line 221 | contains
221      real(kind=dp) :: cutValue
222      
223      scID = SCList%atidToSCType(atomID)
224 <    cutValue = SCList%SCTypes(scID)%sc_rcut
224 >    cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
225    end function getSCCut
226  
227  
# Line 242 | Line 243 | contains
243      if (.not. allocated(MixingMap)) then
244         allocate(MixingMap(nSCtypes, nSCtypes))
245      endif
246 <
246 >    useGeometricDistanceMixing = usesGeometricDistanceMixing()
247      do i = 1, nSCtypes
248  
249         e1 = SCList%SCtypes(i)%epsilon
# Line 271 | Line 272 | contains
272            MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
273            MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* &
274                 (MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
275 +
276            if (i.ne.j) then
277               MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
278               MixingMap(j,i)%m          = MixingMap(i,j)%m
# Line 289 | Line 291 | contains
291  
292    !! routine checks to see if array is allocated, deallocates array if allocated
293    !! and then creates the array to the required size
294 <  subroutine allocateSC(status)
295 <    integer, intent(out) :: status
294 >  subroutine allocateSC()
295 >    integer :: status
296  
297   #ifdef IS_MPI
298      integer :: nAtomsInRow
299      integer :: nAtomsInCol
300   #endif
301      integer :: alloc_stat
300 –
302  
303 +    
304      status = 0
305   #ifdef IS_MPI
306      nAtomsInRow = getNatomsInRow(plan_atom_row)
# Line 311 | Line 313 | contains
313      allocate(frho(nlocal),stat=alloc_stat)
314      if (alloc_stat /= 0) then
315         status = -1
314 –       return
316      end if
317  
318      if (allocated(rho)) deallocate(rho)
319      allocate(rho(nlocal),stat=alloc_stat)
320      if (alloc_stat /= 0) then
321         status = -1
321 –       return
322      end if
323  
324      if (allocated(dfrhodrho)) deallocate(dfrhodrho)
325      allocate(dfrhodrho(nlocal),stat=alloc_stat)
326      if (alloc_stat /= 0) then
327         status = -1
328 –       return
328      end if
329  
330   #ifdef IS_MPI
# Line 334 | Line 333 | contains
333      allocate(rho_tmp(nlocal),stat=alloc_stat)
334      if (alloc_stat /= 0) then
335         status = -1
337 –       return
336      end if
337  
338  
# Line 342 | Line 340 | contains
340      allocate(frho_row(nAtomsInRow),stat=alloc_stat)
341      if (alloc_stat /= 0) then
342         status = -1
345 –       return
343      end if
344      if (allocated(rho_row)) deallocate(rho_row)
345      allocate(rho_row(nAtomsInRow),stat=alloc_stat)
346      if (alloc_stat /= 0) then
347         status = -1
351 –       return
348      end if
349      if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
350      allocate(dfrhodrho_row(nAtomsInRow),stat=alloc_stat)
351      if (alloc_stat /= 0) then
352         status = -1
357 –       return
353      end if
354  
355  
# Line 364 | Line 359 | contains
359      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
360      if (alloc_stat /= 0) then
361         status = -1
367 –       return
362      end if
363      if (allocated(rho_col)) deallocate(rho_col)
364      allocate(rho_col(nAtomsInCol),stat=alloc_stat)
365      if (alloc_stat /= 0) then
366         status = -1
373 –       return
367      end if
368      if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
369      allocate(dfrhodrho_col(nAtomsInCol),stat=alloc_stat)
370      if (alloc_stat /= 0) then
371         status = -1
379 –       return
372      end if
373  
374   #endif
375 <
375 >    if (status == -1) then
376 >       call handleError("SuttonChen:allocateSC","Error in allocating SC arrays")
377 >    end if
378 >    arraysAllocated = .true.
379    end subroutine allocateSC
380  
381    !! C sets rcut to be the largest cutoff of any atype
# Line 394 | Line 389 | contains
389      SC_rcut = rcut
390  
391    end subroutine setCutoffSC
397 –
398 –  subroutine useGeometricMixing()
399 –    useGeometricDistanceMixing = .true.
400 –    haveMixingMap = .false.
401 –    return
402 –  end subroutine useGeometricMixing
403 –  
404 –
392  
393 <
407 <
408 <
409 <
410 <
411 <
393 > !! This array allocates module arrays if needed and builds mixing map.
394    subroutine clean_SC()
395 <
395 >    if (.not.arraysAllocated) call allocateSC()
396 >    if (.not.haveMixingMap) call createMixingMap()
397      ! clean non-IS_MPI first
398      frho = 0.0_dp
399      rho  = 0.0_dp
# Line 430 | Line 413 | contains
413  
414  
415    !! Calculates rho_r
416 <  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq)
416 >  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq, rcut)
417      integer :: atom1,atom2
418      real(kind = dp), dimension(3) :: d
419      real(kind = dp), intent(inout)               :: r
420 <    real(kind = dp), intent(inout)               :: rijsq
420 >    real(kind = dp), intent(inout)               :: rijsq, rcut
421      ! value of electron density rho do to atom i at atom j
422      real(kind = dp) :: rho_i_at_j
423      ! value of electron density rho do to atom j at atom i
# Line 451 | Line 434 | contains
434  
435      ! check to see if we need to be cleaned at the start of a force loop
436  
437 +    if (cleanArrays) call clean_SC()
438 +    cleanArrays = .false.
439  
455 –
456 –
440   #ifdef IS_MPI
441      Atid1 = Atid_row(Atom1)
442      Atid2 = Atid_col(Atom2)
# Line 488 | Line 471 | contains
471      real(kind=dp) :: pot
472      integer :: i,j
473      integer :: atom
491 –    real(kind=dp) :: U,U1,U2
474      integer :: atype1
475      integer :: atid1
476      integer :: myid
477  
478  
497 –    cleanme = .true.
479      !! Scatter the electron density from  pre-pair calculation back to local atoms
480   #ifdef IS_MPI
481      call scatter(rho_row,rho,plan_atom_row,sc_err)
# Line 518 | Line 499 | contains
499      do atom = 1, nlocal
500         Myid = SCList%atidtoSctype(Atid(atom))
501         frho(atom) = -SCList%SCTypes(Myid)%c * &
502 <            SCList%SCTypes(Myid)%epsilon * sqrt(rho(i))
502 >            SCList%SCTypes(Myid)%epsilon * sqrt(rho(atom))
503  
504         dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
505 <       pot = pot + u
505 >
506 >       pot = pot + frho(atom)
507      enddo
508  
509   #ifdef IS_MPI
# Line 549 | Line 531 | contains
531  
532  
533    !! Does Sutton-Chen  pairwise Force calculation.  
534 <  subroutine do_sc_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
534 >  subroutine do_sc_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
535         pot, f, do_pot)
536      !Arguments    
537      integer, intent(in) ::  atom1, atom2
538 <    real( kind = dp ), intent(in) :: rij, r2
538 >    real( kind = dp ), intent(in) :: rij, r2, rcut
539      real( kind = dp ) :: pot, sw, vpair
540      real( kind = dp ), dimension(3,nLocal) :: f
541      real( kind = dp ), intent(in), dimension(3) :: d
# Line 565 | Line 547 | contains
547      real( kind = dp ) :: dvpdr
548      real( kind = dp ) :: drhodr
549      real( kind = dp ) :: dudr
568 –    real( kind = dp ) :: rcij
550      real( kind = dp ) :: drhoidr,drhojdr
551      real( kind = dp ) :: Fx,Fy,Fz
552 <    real( kind = dp ) :: r,d2pha,phb,d2phb
552 >    real( kind = dp ) :: d2pha,phb,d2phb
553      real( kind = dp ) :: pot_temp,vptmp
554      real( kind = dp ) :: epsilonij,aij,nij,mij,vcij
555      integer :: id1,id2
# Line 578 | Line 559 | contains
559      !Local Variables
560  
561      ! write(*,*) "Frho: ", Frho(atom1)
562 +    
563 +    cleanArrays = .true.
564  
582 –
565      dvpdr = 0.0E0_DP
566  
567  
# Line 605 | Line 587 | contains
587         mij       = MixingMap(mytype_atom1,mytype_atom2)%m
588         vcij      = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot              
589  
590 <       vptmp = epsilonij*((aij/r)**nij)
590 >       vptmp = epsilonij*((aij/rij)**nij)
591  
592  
593 <       dvpdr = -nij*vptmp/r
594 <       drhodr = -mij*((aij/r)**mij)/r
593 >       dvpdr = -nij*vptmp/rij
594 >       drhodr = -mij*((aij/rij)**mij)/rij
595  
596        
597         dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &

Diff Legend

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