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 2530 by chuckv, Fri Dec 30 00:18:28 2005 UTC vs.
Revision 2717 by gezelter, Mon Apr 17 21:49:12 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 50 | Line 50 | module suttonchen
50    use atype_module
51    use vector_class
52    use fForceOptions
53 +  use interpolation
54   #ifdef IS_MPI
55    use mpiSimulation
56   #endif
# Line 59 | Line 60 | module suttonchen
60   #include "UseTheForce/DarkSide/fInteractionMap.h"
61  
62    INTEGER, PARAMETER :: DP = selected_real_kind(15)
63 +  !! number of points for the spline approximations
64 +  INTEGER, PARAMETER :: np = 3000
65  
66    logical, save :: SC_FF_initialized = .false.
67    integer, save :: SC_Mixing_Policy
# Line 66 | Line 69 | module suttonchen
69    logical, save :: haveRcut = .false.
70    logical, save :: haveMixingMap = .false.
71    logical, save :: useGeometricDistanceMixing = .false.
72 +  logical, save :: cleanArrays = .true.
73 +  logical, save :: arraysAllocated = .false.
74  
75  
71 –
72 –
76    character(len = statusMsgSize) :: errMesg
77    integer :: sc_err
78  
79    character(len = 200) :: errMsg
80    character(len=*), parameter :: RoutineName =  "Sutton-Chen MODULE"
81 <  !! Logical that determines if eam arrays should be zeroed
79 <  logical :: cleanme = .true.
80 <  logical :: nmflag  = .false.
81 <
82 <
81 >  
82    type, private :: SCtype
83 <     integer           :: atid      
84 <     real(kind=dp)     :: c
85 <     real(kind=dp)     :: m
86 <     real(kind=dp)     :: n
87 <     real(kind=dp)     :: alpha
88 <     real(kind=dp)     :: epsilon
89 <     real(kind=dp)     :: sc_rcut
83 >     integer       :: atid      
84 >     real(kind=dp) :: c
85 >     real(kind=dp) :: m
86 >     real(kind=dp) :: n
87 >     real(kind=dp) :: alpha
88 >     real(kind=dp) :: epsilon
89 >     real(kind=dp) :: sc_rcut
90    end type SCtype
91 +  
92  
93 –
93    !! Arrays for derivatives used in force calculation
94    real( kind = dp), dimension(:), allocatable :: rho
95    real( kind = dp), dimension(:), allocatable :: frho
96    real( kind = dp), dimension(:), allocatable :: dfrhodrho
97 <
99 <
100 <
97 >  
98    !! Arrays for MPI storage
99   #ifdef IS_MPI
100    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
# Line 111 | Line 108 | module suttonchen
108  
109    type, private :: SCTypeList
110       integer           :: nSCTypes = 0
111 <     integer           :: currentSCtype = 0
112 <
113 <     type (SCtype), pointer  :: SCtypes(:) => null()
117 <     integer, pointer         :: atidToSCtype(:) => null()
111 >     integer           :: currentSCtype = 0    
112 >     type (SCtype), pointer :: SCtypes(:) => null()
113 >     integer, pointer       :: atidToSCtype(:) => null()
114    end type SCTypeList
115  
116    type (SCTypeList), save :: SCList
117  
118 <
123 <
124 <
125 < type :: MixParameters
118 >  type:: MixParameters
119       real(kind=DP) :: alpha
120       real(kind=DP) :: epsilon
121 <     real(kind=DP) :: m
121 >     real(kind=DP) :: m
122       real(Kind=DP) :: n
130 –     real(kind=DP) :: vpair_pot
123       real(kind=dp) :: rCut
124 +     real(kind=dp) :: vCut
125       logical       :: rCutWasSet = .false.
126 <
126 >     type(cubicSpline) :: V
127 >     type(cubicSpline) :: phi
128    end type MixParameters
129  
130    type(MixParameters), dimension(:,:), allocatable :: MixingMap
131  
138 –
139 –
132    public :: setCutoffSC
133    public :: do_SC_pair
134    public :: newSCtype
# Line 158 | Line 150 | contains
150      real (kind = dp )                      :: n ! Pair Potential Exponent
151      real (kind = dp )                      :: alpha ! Length Scaling
152      real (kind = dp )                      :: epsilon ! Energy Scaling
161 –
162 –
153      integer                                :: c_ident
154      integer                                :: status
165 –
155      integer                                :: nAtypes,nSCTypes,myATID
156      integer                                :: maxVals
157      integer                                :: alloc_stat
# Line 189 | Line 178 | contains
178  
179      myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
180      SCList%atidToSCType(myATID) = current
192 –
193 –
181    
182      SCList%SCTypes(current)%atid         = c_ident
183      SCList%SCTypes(current)%alpha        = alpha
# Line 210 | Line 197 | contains
197         deallocate(SCList%atidToSCtype)
198         SCList%atidToSCtype=>null()
199      end if
200 + ! Reset Capacity
201 +    SCList%nSCTypes = 0
202 +    SCList%currentSCtype=0
203  
214 –
204    end subroutine destroySCTypes
205  
217 –
218 –
206    function getSCCut(atomID) result(cutValue)
207      integer, intent(in) :: atomID
208      integer :: scID
# Line 225 | Line 212 | contains
212      cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
213    end function getSCCut
214  
228 –
229 –
230 –
215    subroutine createMixingMap()
216 <    integer :: nSCtypes, i, j
217 <    real ( kind = dp ) :: e1, e2,m1,m2,alpha1,alpha2,n1,n2
218 <    real ( kind = dp ) :: rcut6, tp6, tp12
219 <    logical :: isSoftCore1, isSoftCore2, doShift
216 >    integer :: nSCtypes, i, j, k
217 >    real ( kind = dp ) :: e1, e2, m1, m2, alpha1, alpha2, n1, n2
218 >    real ( kind = dp ) :: epsilon, m, n, alpha, rCut, vCut, dr, r
219 >    real ( kind = dp ), dimension(np) :: rvals, vvals, phivals
220  
221      if (SCList%currentSCtype == 0) then
222         call handleError("SuttonChen", "No members in SCMap")
# Line 252 | Line 236 | contains
236         n1 = SCList%SCtypes(i)%n
237         alpha1 = SCList%SCtypes(i)%alpha
238  
239 <       do j = i, nSCtypes
239 >       do j = 1, nSCtypes
240            
257 –
241            e2 = SCList%SCtypes(j)%epsilon
242            m2 = SCList%SCtypes(j)%m
243            n2 = SCList%SCtypes(j)%n
244            alpha2 = SCList%SCtypes(j)%alpha
245  
246            if (useGeometricDistanceMixing) then
247 <             MixingMap(i,j)%alpha = sqrt(alpha1 * alpha2) !SC formulation
247 >             alpha = sqrt(alpha1 * alpha2) !SC formulation
248            else
249 <             MixingMap(i,j)%alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
249 >             alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
250            endif
251 +          rcut = 2.0_dp * alpha
252 +          epsilon = sqrt(e1 * e2)
253 +          m = 0.5_dp*(m1+m2)
254 +          n = 0.5_dp*(n1+n2)
255  
256 <          MixingMap(i,j)%epsilon = sqrt(e1 * e2)
257 <          MixingMap(i,j)%m = 0.5_dp*(m1+m2)
258 <          MixingMap(i,j)%n = 0.5_dp*(n1+n2)
259 <          MixingMap(i,j)%alpha = 0.5_dp*(alpha1+alpha2)
260 <          MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
261 <          MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* &
262 <               (MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
263 <          if (i.ne.j) then
264 <             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
265 <             MixingMap(j,i)%m          = MixingMap(i,j)%m
266 <             MixingMap(j,i)%n          = MixingMap(i,j)%n
267 <             MixingMap(j,i)%alpha      = MixingMap(i,j)%alpha
268 <             MixingMap(j,i)%rcut = MixingMap(i,j)%rcut
269 <             MixingMap(j,i)%vpair_pot = MixingMap(i,j)%vpair_pot
270 <          endif
256 >          dr = (rCut) / dble(np-1)
257 >          rvals(1) = 0.0d0
258 >          vvals(1) = 0.0d0
259 >          phivals(1) = 0.0d0
260 >
261 >          do k = 2, np
262 >             r = dble(k-1)*dr
263 >             rvals(k) = r
264 >             vvals(k) = epsilon*((alpha/r)**n)
265 >             phivals(k) = (alpha/r)**m
266 >          enddo
267 >
268 >          vCut = epsilon*((alpha/rCut)**n)
269 >
270 >          call newSpline(MixingMap(i,j)%V, rvals, vvals, &
271 >               0.0d0, 0.0d0, .true.)
272 >          call newSpline(MixingMap(i,j)%phi, rvals, phivals, &
273 >               0.0d0, 0.0d0, .true.)
274 >
275 >          MixingMap(i,j)%epsilon = epsilon
276 >          MixingMap(i,j)%m = m
277 >          MixingMap(i,j)%n = n
278 >          MixingMap(i,j)%alpha = alpha
279 >          MixingMap(i,j)%rCut = rcut
280 >          MixingMap(i,j)%vCut = vCut
281         enddo
282      enddo
283      
# Line 291 | Line 288 | contains
288  
289    !! routine checks to see if array is allocated, deallocates array if allocated
290    !! and then creates the array to the required size
291 <  subroutine allocateSC(status)
292 <    integer, intent(out) :: status
291 >  subroutine allocateSC()
292 >    integer :: status
293  
294   #ifdef IS_MPI
295      integer :: nAtomsInRow
# Line 300 | Line 297 | contains
297   #endif
298      integer :: alloc_stat
299  
300 <
300 >    
301      status = 0
302   #ifdef IS_MPI
303      nAtomsInRow = getNatomsInRow(plan_atom_row)
304      nAtomsInCol = getNatomsInCol(plan_atom_col)
305   #endif
306  
310 –
311 –
307      if (allocated(frho)) deallocate(frho)
308      allocate(frho(nlocal),stat=alloc_stat)
309      if (alloc_stat /= 0) then
310         status = -1
316 –       return
311      end if
312  
313      if (allocated(rho)) deallocate(rho)
314      allocate(rho(nlocal),stat=alloc_stat)
315      if (alloc_stat /= 0) then
316         status = -1
323 –       return
317      end if
318  
319      if (allocated(dfrhodrho)) deallocate(dfrhodrho)
320      allocate(dfrhodrho(nlocal),stat=alloc_stat)
321      if (alloc_stat /= 0) then
322         status = -1
330 –       return
323      end if
324  
325   #ifdef IS_MPI
# Line 336 | Line 328 | contains
328      allocate(rho_tmp(nlocal),stat=alloc_stat)
329      if (alloc_stat /= 0) then
330         status = -1
339 –       return
331      end if
332  
333  
# Line 344 | Line 335 | contains
335      allocate(frho_row(nAtomsInRow),stat=alloc_stat)
336      if (alloc_stat /= 0) then
337         status = -1
347 –       return
338      end if
339      if (allocated(rho_row)) deallocate(rho_row)
340      allocate(rho_row(nAtomsInRow),stat=alloc_stat)
341      if (alloc_stat /= 0) then
342         status = -1
353 –       return
343      end if
344      if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
345      allocate(dfrhodrho_row(nAtomsInRow),stat=alloc_stat)
346      if (alloc_stat /= 0) then
347         status = -1
359 –       return
348      end if
349  
350  
# Line 366 | Line 354 | contains
354      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
355      if (alloc_stat /= 0) then
356         status = -1
369 –       return
357      end if
358      if (allocated(rho_col)) deallocate(rho_col)
359      allocate(rho_col(nAtomsInCol),stat=alloc_stat)
360      if (alloc_stat /= 0) then
361         status = -1
375 –       return
362      end if
363      if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
364      allocate(dfrhodrho_col(nAtomsInCol),stat=alloc_stat)
365      if (alloc_stat /= 0) then
366         status = -1
381 –       return
367      end if
368  
369   #endif
370 <
370 >    if (status == -1) then
371 >       call handleError("SuttonChen:allocateSC","Error in allocating SC arrays")
372 >    end if
373 >    arraysAllocated = .true.
374    end subroutine allocateSC
375  
376 <  !! C sets rcut to be the largest cutoff of any atype
377 <  !! present in this simulation. Doesn't include all atypes
390 <  !! sim knows about, just those in the simulation.
391 <  subroutine setCutoffSC(rcut, status)
392 <    real(kind=dp) :: rcut
393 <    integer :: status
394 <    status = 0
395 <
376 >  subroutine setCutoffSC(rcut)
377 >    real(kind=dp) :: rcut
378      SC_rcut = rcut
397 –
379    end subroutine setCutoffSC
380 <
380 >  
381 >  !! This array allocates module arrays if needed and builds mixing map.
382    subroutine clean_SC()
383 <
383 >    if (.not.arraysAllocated) call allocateSC()
384 >    if (.not.haveMixingMap) call createMixingMap()
385      ! clean non-IS_MPI first
386      frho = 0.0_dp
387      rho  = 0.0_dp
# Line 415 | Line 398 | contains
398   #endif
399    end subroutine clean_SC
400  
418 –
419 –
401    !! Calculates rho_r
402 <  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq, rcut)
402 >  subroutine calc_sc_prepair_rho(atom1, atom2, d, r, rijsq, rcut)
403      integer :: atom1,atom2
404      real(kind = dp), dimension(3) :: d
405      real(kind = dp), intent(inout)               :: r
# Line 439 | Line 420 | contains
420  
421      ! check to see if we need to be cleaned at the start of a force loop
422  
423 +    if (cleanArrays) call clean_SC()
424 +    cleanArrays = .false.
425  
443 –
444 –
426   #ifdef IS_MPI
427      Atid1 = Atid_row(Atom1)
428      Atid2 = Atid_col(Atom2)
# Line 452 | Line 433 | contains
433  
434      Myid_atom1 = SCList%atidtoSCtype(Atid1)
435      Myid_atom2 = SCList%atidtoSCtype(Atid2)
436 +    
437 +    call lookupUniformSpline(MixingMap(myid_atom1,myid_atom2)%phi, r, &
438 +         rho_i_at_j)
439 +    rho_j_at_i = rho_i_at_j
440  
456 –
457 –
458 –       rho_i_at_j = (MixingMap(Myid_atom1,Myid_atom2)%alpha/r)&
459 –            **MixingMap(Myid_atom1,Myid_atom2)%m
460 –       rho_j_at_i = rho_i_at_j
461 –
441   #ifdef  IS_MPI
442 <       rho_col(atom2) = rho_col(atom2) + rho_i_at_j
443 <       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
442 >    rho_col(atom2) = rho_col(atom2) + rho_i_at_j
443 >    rho_row(atom1) = rho_row(atom1) + rho_j_at_i
444   #else
445 <       rho(atom2) = rho(atom2) + rho_i_at_j
446 <       rho(atom1) = rho(atom1) + rho_j_at_i
445 >    rho(atom2) = rho(atom2) + rho_i_at_j
446 >    rho(atom1) = rho(atom1) + rho_j_at_i
447   #endif
448 <
448 >    
449    end subroutine calc_sc_prepair_rho
450  
451  
452 < !! Calculate the rho_a for all local atoms
452 >  !! Calculate the rho_a for all local atoms
453    subroutine calc_sc_preforce_Frho(nlocal,pot)
454      integer :: nlocal
455      real(kind=dp) :: pot
456      integer :: i,j
457      integer :: atom
479 –    real(kind=dp) :: U,U1,U2
458      integer :: atype1
459      integer :: atid1
460      integer :: myid
461  
462 <
463 <    cleanme = .true.
486 <    !! Scatter the electron density from  pre-pair calculation back to local atoms
462 >    !! Scatter the electron density from  pre-pair calculation back to
463 >    !! local atoms
464   #ifdef IS_MPI
465      call scatter(rho_row,rho,plan_atom_row,sc_err)
466      if (sc_err /= 0 ) then
# Line 499 | Line 476 | contains
476  
477      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
478   #endif
479 <
503 <
504 <
479 >    
480      !! Calculate F(rho) and derivative
481      do atom = 1, nlocal
482         Myid = SCList%atidtoSctype(Atid(atom))
483 <       frho(atom) = -SCList%SCTypes(Myid)%c * &
484 <            SCList%SCTypes(Myid)%epsilon * sqrt(rho(i))
483 >       frho(atom) = - SCList%SCTypes(Myid)%c * &
484 >            SCList%SCTypes(Myid)%epsilon * sqrt(rho(atom))
485  
486         dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
487 <       pot = pot + u
487 >
488 >       pot = pot + frho(atom)
489      enddo
490  
491   #ifdef IS_MPI
# Line 533 | Line 509 | contains
509   #endif
510      
511      
512 <  end subroutine calc_sc_preforce_Frho
513 <
538 <
512 >  end subroutine calc_sc_preforce_Frho  
513 >  
514    !! Does Sutton-Chen  pairwise Force calculation.  
515    subroutine do_sc_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
516         pot, f, do_pot)
# Line 549 | Line 524 | contains
524  
525      logical, intent(in) :: do_pot
526  
527 <    real( kind = dp ) :: drdx,drdy,drdz
527 >    real( kind = dp ) :: drdx, drdy, drdz
528      real( kind = dp ) :: dvpdr
529 <    real( kind = dp ) :: drhodr
529 >    real( kind = dp ) :: rhtmp, drhodr
530      real( kind = dp ) :: dudr
556 –    real( kind = dp ) :: rcij
557 –    real( kind = dp ) :: drhoidr,drhojdr
531      real( kind = dp ) :: Fx,Fy,Fz
532 <    real( kind = dp ) :: r,d2pha,phb,d2phb
533 <    real( kind = dp ) :: pot_temp,vptmp
534 <    real( kind = dp ) :: epsilonij,aij,nij,mij,vcij
562 <    integer :: id1,id2
532 >    real( kind = dp ) :: pot_temp, vptmp
533 >    real( kind = dp ) :: rcij, vcij
534 >    integer :: id1, id2
535      integer  :: mytype_atom1
536      integer  :: mytype_atom2
537 <    integer  :: atid1,atid2
537 >    integer  :: atid1, atid2
538      !Local Variables
539 +    
540 +    cleanArrays = .true.
541  
568 –    ! write(*,*) "Frho: ", Frho(atom1)
569 –
570 –
571 –    dvpdr = 0.0E0_DP
572 –
573 –
542   #ifdef IS_MPI
543         atid1 = atid_row(atom1)
544         atid2 = atid_col(atom2)
# Line 585 | Line 553 | contains
553         drdx = d(1)/rij
554         drdy = d(2)/rij
555         drdz = d(3)/rij
588 –
556                  
557 <       epsilonij = MixingMap(mytype_atom1,mytype_atom2)%epsilon
558 <       aij       = MixingMap(mytype_atom1,mytype_atom2)%alpha
592 <       nij       = MixingMap(mytype_atom1,mytype_atom2)%n
593 <       mij       = MixingMap(mytype_atom1,mytype_atom2)%m
594 <       vcij      = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot              
595 <
596 <       vptmp = epsilonij*((aij/r)**nij)
597 <
598 <
599 <       dvpdr = -nij*vptmp/r
600 <       drhodr = -mij*((aij/r)**mij)/r
601 <
557 >       rcij = MixingMap(mytype_atom1,mytype_atom2)%rCut
558 >       vcij = MixingMap(mytype_atom1,mytype_atom2)%vCut
559        
560 <       dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
561 <            + dvpdr
605 <      
606 <       pot_temp = vptmp + vcij
607 <
560 >       call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%phi, &
561 >            rij, rhtmp, drhodr)
562  
563 +       call lookupUniformSpline1d(MixingMap(mytype_atom1, mytype_atom2)%V, &
564 +            rij, vptmp, dvpdr)
565 +      
566   #ifdef IS_MPI
567 <       dudr = drhodr*(dfrhodrho_row(atom1)+dfrhodrho_col(atom2)) &
611 <            + dvpdr
612 <
567 >       dudr = drhodr*(dfrhodrho_row(atom1) + dfrhodrho_col(atom2)) + dvpdr
568   #else
569 <       dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
615 <            + dvpdr
569 >       dudr = drhodr*(dfrhodrho(atom1)+ dfrhodrho(atom2)) + dvpdr
570   #endif
571 <
572 <
571 >              
572 >       pot_temp = vptmp - vcij
573 >
574         fx = dudr * drdx
575         fy = dudr * drdy
576         fz = dudr * drdz
577  
623 –
578   #ifdef IS_MPI
579         if (do_pot) then
580            pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + (pot_temp)*0.5
# Line 665 | Line 619 | contains
619            fpair(3) = fpair(3) + fz
620  
621         endif
668 –
669 –
622    end subroutine do_sc_pair
671 –
672 –
673 –
623   end module suttonchen

Diff Legend

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