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 2680 by chuckv, Thu Mar 30 21:08:48 2006 UTC vs.
Revision 2717 by gezelter, Mon Apr 17 21:49:12 2006 UTC

# 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 75 | Line 78 | module suttonchen
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 :: nmflag  = .false.
80 <
81 <
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  
92 –
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 <
98 <
99 <
97 >  
98    !! Arrays for MPI storage
99   #ifdef IS_MPI
100    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
# Line 110 | Line 108 | module suttonchen
108  
109    type, private :: SCTypeList
110       integer           :: nSCTypes = 0
111 <     integer           :: currentSCtype = 0
112 <
113 <     type (SCtype), pointer  :: SCtypes(:) => null()
116 <     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 <
122 <
123 <
124 < 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
129 –     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  
137 –
138 –
132    public :: setCutoffSC
133    public :: do_SC_pair
134    public :: newSCtype
# Line 157 | 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
160 –
161 –
153      integer                                :: c_ident
154      integer                                :: status
164 –
155      integer                                :: nAtypes,nSCTypes,myATID
156      integer                                :: maxVals
157      integer                                :: alloc_stat
# Line 188 | Line 178 | contains
178  
179      myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
180      SCList%atidToSCType(myATID) = current
191 –
192 –
181    
182      SCList%SCTypes(current)%atid         = c_ident
183      SCList%SCTypes(current)%alpha        = alpha
# Line 210 | Line 198 | contains
198         SCList%atidToSCtype=>null()
199      end if
200   ! Reset Capacity
201 <    SCList% nSCTypes = 0
201 >    SCList%nSCTypes = 0
202      SCList%currentSCtype=0
203  
204    end subroutine destroySCTypes
205  
218 –
219 –
206    function getSCCut(atomID) result(cutValue)
207      integer, intent(in) :: atomID
208      integer :: scID
# Line 226 | Line 212 | contains
212      cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
213    end function getSCCut
214  
229 –
230 –
231 –
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 253 | 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            
258 –
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)
274 <          MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
275 <          MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* &
276 <               (MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
256 >          dr = (rCut) / dble(np-1)
257 >          rvals(1) = 0.0d0
258 >          vvals(1) = 0.0d0
259 >          phivals(1) = 0.0d0
260  
261 <          if (i.ne.j) then
262 <             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
263 <             MixingMap(j,i)%m          = MixingMap(i,j)%m
264 <             MixingMap(j,i)%n          = MixingMap(i,j)%n
265 <             MixingMap(j,i)%alpha      = MixingMap(i,j)%alpha
266 <             MixingMap(j,i)%rcut = MixingMap(i,j)%rcut
267 <             MixingMap(j,i)%vpair_pot = MixingMap(i,j)%vpair_pot
268 <          endif
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 309 | Line 304 | contains
304      nAtomsInCol = getNatomsInCol(plan_atom_col)
305   #endif
306  
312 –
313 –
307      if (allocated(frho)) deallocate(frho)
308      allocate(frho(nlocal),stat=alloc_stat)
309      if (alloc_stat /= 0) then
# Line 380 | Line 373 | contains
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
385 <  !! sim knows about, just those in the simulation.
386 <  subroutine setCutoffSC(rcut, status)
387 <    real(kind=dp) :: rcut
388 <    integer :: status
389 <    status = 0
390 <
376 >  subroutine setCutoffSC(rcut)
377 >    real(kind=dp) :: rcut
378      SC_rcut = rcut
392 –
379    end subroutine setCutoffSC
380 <
381 < !! This array allocates module arrays if needed and builds mixing map.
380 >  
381 >  !! This array allocates module arrays if needed and builds mixing map.
382    subroutine clean_SC()
383      if (.not.arraysAllocated) call allocateSC()
384      if (.not.haveMixingMap) call createMixingMap()
# Line 412 | Line 398 | contains
398   #endif
399    end subroutine clean_SC
400  
415 –
416 –
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 449 | 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  
453 –
454 –
455 –       rho_i_at_j = (MixingMap(Myid_atom1,Myid_atom2)%alpha/r)&
456 –            **MixingMap(Myid_atom1,Myid_atom2)%m
457 –       rho_j_at_i = rho_i_at_j
458 –
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
# Line 477 | Line 459 | contains
459      integer :: atid1
460      integer :: myid
461  
462 <
463 <    !! 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 494 | Line 476 | contains
476  
477      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
478   #endif
479 <
498 <
499 <
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 * &
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)
# Line 529 | Line 509 | contains
509   #endif
510      
511      
512 <  end subroutine calc_sc_preforce_Frho
513 <
534 <
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 545 | 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
552 –    real( kind = dp ) :: drhoidr,drhojdr
531      real( kind = dp ) :: Fx,Fy,Fz
532 <    real( kind = dp ) :: d2pha,phb,d2phb
533 <    real( kind = dp ) :: pot_temp,vptmp
534 <    real( kind = dp ) :: epsilonij,aij,nij,mij,vcij
557 <    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
562 –
563 –    ! write(*,*) "Frho: ", Frho(atom1)
539      
540      cleanArrays = .true.
541  
567 –    dvpdr = 0.0E0_DP
568 –
569 –
542   #ifdef IS_MPI
543         atid1 = atid_row(atom1)
544         atid2 = atid_col(atom2)
# Line 581 | Line 553 | contains
553         drdx = d(1)/rij
554         drdy = d(2)/rij
555         drdz = d(3)/rij
584 –
556                  
557 <       epsilonij = MixingMap(mytype_atom1,mytype_atom2)%epsilon
558 <       aij       = MixingMap(mytype_atom1,mytype_atom2)%alpha
588 <       nij       = MixingMap(mytype_atom1,mytype_atom2)%n
589 <       mij       = MixingMap(mytype_atom1,mytype_atom2)%m
590 <       vcij      = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot              
591 <
592 <       vptmp = epsilonij*((aij/rij)**nij)
593 <
594 <
595 <       dvpdr = -nij*vptmp/rij
596 <       drhodr = -mij*((aij/rij)**mij)/rij
597 <
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
601 <      
602 <       pot_temp = vptmp + vcij
603 <
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)) &
607 <            + dvpdr
608 <
567 >       dudr = drhodr*(dfrhodrho_row(atom1) + dfrhodrho_col(atom2)) + dvpdr
568   #else
569 <       dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
611 <            + 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  
619 –
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 661 | Line 619 | contains
619            fpair(3) = fpair(3) + fz
620  
621         endif
664 –
665 –
622    end subroutine do_sc_pair
667 –
668 –
669 –
623   end module suttonchen

Diff Legend

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