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

Comparing trunk/OOPSE-3.0/src/UseTheForce/DarkSide/suttonchen.F90 (file contents):
Revision 2426 by chuckv, Thu Nov 3 23:06:00 2005 UTC vs.
Revision 2427 by chuckv, Mon Nov 14 21:29:36 2005 UTC

# Line 63 | Line 63 | module suttonchen
63    integer, save :: SC_Mixing_Policy
64    real(kind = dp), save :: SC_rcut
65    logical, save :: haveRcut = .false.
66 <  logical, save,:: haveMixingMap = .false.
66 >  logical, save :: haveMixingMap = .false.
67 >  logical, save :: useGeometricDistanceMixing = .false.
68  
69 +
70 +
71 +
72    character(len = statusMsgSize) :: errMesg
73    integer :: sc_err
74  
# Line 82 | Line 86 | module suttonchen
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  
# Line 89 | Line 94 | module suttonchen
94    real( kind = dp), dimension(:), allocatable :: rho
95    real( kind = dp), dimension(:), allocatable :: frho
96    real( kind = dp), dimension(:), allocatable :: dfrhodrho
92 –  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
97  
98  
99 +
100    !! Arrays for MPI storage
101   #ifdef IS_MPI
102    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
# Line 101 | Line 106 | module suttonchen
106    real( kind = dp),save, dimension(:), allocatable :: rho_row
107    real( kind = dp),save, dimension(:), allocatable :: rho_col
108    real( kind = dp),save, dimension(:), allocatable :: rho_tmp
104 –  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col
105 –  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
109   #endif
110  
111    type, private :: SCTypeList
# Line 121 | Line 124 | module suttonchen
124   type :: MixParameters
125       real(kind=DP) :: alpha
126       real(kind=DP) :: epsilon
127 <     real(kind=dp) :: sigma6
127 >     real(kind=DP) :: m
128 >     real(Kind=DP) :: n
129 >     real(kind=DP) :: vpair_pot
130       real(kind=dp) :: rCut
126 –     real(kind=dp) :: vpair_pot
131       logical       :: rCutWasSet = .false.
132 <     logical       :: shiftedPot
129 <     logical       :: isSoftCore = .false.
132 >
133    end type MixParameters
134  
135    type(MixParameters), dimension(:,:), allocatable :: MixingMap
136  
137  
138  
136 –  public :: init_SC_FF
139    public :: setCutoffSC
140    public :: do_SC_pair
141    public :: newSCtype
# Line 141 | Line 143 | module suttonchen
143    public :: clean_SC
144    public :: destroySCtypes
145    public :: getSCCut
146 <  public :: setSCDefaultCutoff
147 <  public :: setSCUniformCutoff
146 > ! public :: setSCDefaultCutoff
147 > ! public :: setSCUniformCutoff
148 >  public :: useGeometricMixing
149  
150   contains
151  
152  
153 <  subroutine newSCtype(c,m,n,alpha,epsilon,c_ident,status)
154 <    real (kind = dp )                      :: lattice_constant
155 <    integer                                :: eam_nrho
156 <    real (kind = dp )                      :: eam_drho
157 <    integer                                :: eam_nr
158 <    real (kind = dp )                      :: eam_dr
159 <    real (kind = dp )                      :: rcut
160 <    real (kind = dp ), dimension(eam_nr)   :: eam_Z_r
158 <    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
159 <    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
153 >  subroutine newSCtype(c_ident,c,m,n,alpha,epsilon,status)
154 >    real (kind = dp )                      :: c ! Density Scaling
155 >    real (kind = dp )                      :: m ! Density Exponent
156 >    real (kind = dp )                      :: n ! Pair Potential Exponent
157 >    real (kind = dp )                      :: alpha ! Length Scaling
158 >    real (kind = dp )                      :: epsilon ! Energy Scaling
159 >
160 >
161      integer                                :: c_ident
162      integer                                :: status
163  
# Line 170 | Line 171 | contains
171  
172  
173      !! Assume that atypes has already been set and get the total number of types in atypes
173 –    !! Also assume that every member of atypes is a EAM model.
174  
175  
176      ! check to see if this is the first time into
177 <    if (.not.associated(SCTypeList%EAMParams)) then
177 >    if (.not.associated(SCList%SCTypes)) then
178         call getMatchingElementList(atypes, "is_SuttonChen", .true., nSCtypes, MatchList)
179 <       SCTypeList%nSCtypes = nSCtypes
180 <       allocate(SCTypeList%SCTypes(nSCTypes))
179 >       SCList%nSCtypes = nSCtypes
180 >       allocate(SCList%SCTypes(nSCTypes))
181         nAtypes = getSize(atypes)
182 <       allocate(SCTypeList%atidToSCType(nAtypes))
182 >       allocate(SCList%atidToSCType(nAtypes))
183      end if
184  
185 <    SCTypeList%currentSCType = SCTypeList%currentSCType + 1
186 <    current = SCTypeList%currentSCType
185 >    SCList%currentSCType = SCList%currentSCType + 1
186 >    current = SCList%currentSCType
187  
188      myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
189 <    SCTypeList%atidToSCType(myATID) = current
189 >    SCList%atidToSCType(myATID) = current
190  
191  
192    
193 <    SCTypeList%SCTypes(current)%atid         = c_ident
194 <    SCTypeList%SCTypes(current)%alpha        = alpha
195 <    SCTypeList%SCTypes(current)%c            = c
196 <    SCTypeList%SCTypes(current)%m            = m
197 <    SCTypeList%SCTypes(current)%n            = n
198 <    SCTypeList%SCTypes(current)%epsilon      = epsilon
193 >    SCList%SCTypes(current)%atid         = c_ident
194 >    SCList%SCTypes(current)%alpha        = alpha
195 >    SCList%SCTypes(current)%c            = c
196 >    SCList%SCTypes(current)%m            = m
197 >    SCList%SCTypes(current)%n            = n
198 >    SCList%SCTypes(current)%epsilon      = epsilon
199    end subroutine newSCtype
200  
201    
202    subroutine destroySCTypes()
203 +    if (associated(SCList%SCtypes)) then
204 +       deallocate(SCList%SCTypes)
205 +       SCList%SCTypes=>null()
206 +    end if
207 +    if (associated(SCList%atidToSCtype)) then
208 +       deallocate(SCList%atidToSCtype)
209 +       SCList%atidToSCtype=>null()
210 +    end if
211  
204 –    
205 –    if(allocated(SCTypeList)) deallocate(SCTypeList)
212  
207 –
208 –
213    end subroutine destroySCTypes
214  
215  
216  
217    function getSCCut(atomID) result(cutValue)
218      integer, intent(in) :: atomID
219 <    integer :: eamID
219 >    integer :: scID
220      real(kind=dp) :: cutValue
221      
222 <    eamID = SCTypeList%atidToEAMType(atomID)
223 <    cutValue = SCTypeList%EAMParams(eamID)%eam_rcut
222 >    scID = SCList%atidToSCType(atomID)
223 >    cutValue = SCList%SCTypes(scID)%sc_rcut
224    end function getSCCut
225  
226  
# Line 254 | Line 258 | contains
258            n2 = SCList%SCtypes(j)%n
259            alpha2 = SCList%SCtypes(j)%alpha
260  
261 <          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
261 >          if (useGeometricDistanceMixing) then
262 >             MixingMap(i,j)%alpha = sqrt(alpha1 * alpha2) !SC formulation
263 >          else
264 >             MixingMap(i,j)%alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
265 >          endif
266 >
267 >          MixingMap(i,j)%epsilon = sqrt(e1 * e2)
268            MixingMap(i,j)%m = 0.5_dp*(m1+m2)
269            MixingMap(i,j)%n = 0.5_dp*(n1+n2)
270            MixingMap(i,j)%alpha = 0.5_dp*(alpha1+alpha2)
271            MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
272 <          MixingMap(i,j)%vpair_rcut = MixingMap%epsilon(i,j)* &
273 <               (MixingMap%alpha(i,j)/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
274 <
272 >          MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* &
273 >               (MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
274 >          if (i.ne.j) then
275 >             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
276 >             MixingMap(j,i)%m          = MixingMap(i,j)%m
277 >             MixingMap(j,i)%n          = MixingMap(i,j)%n
278 >             MixingMap(j,i)%alpha      = MixingMap(i,j)%alpha
279 >             MixingMap(j,i)%rcut = MixingMap(i,j)%rcut
280 >             MixingMap(j,i)%vpair_pot = MixingMap(i,j)%vpair_pot
281 >          endif
282         enddo
283      enddo
284      
# Line 311 | Line 328 | contains
328         return
329      end if
330  
314 –    if (allocated(d2frhodrhodrho)) deallocate(d2frhodrhodrho)
315 –    allocate(d2frhodrhodrho(nlocal),stat=alloc_stat)
316 –    if (alloc_stat /= 0) then
317 –       status = -1
318 –       return
319 –    end if
320 –
331   #ifdef IS_MPI
332  
333      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 342 | Line 352 | contains
352      end if
353      if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
354      allocate(dfrhodrho_row(nAtomsInRow),stat=alloc_stat)
345 –    if (alloc_stat /= 0) then
346 –       status = -1
347 –       return
348 –    end if
349 –    if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
350 –    allocate(d2frhodrhodrho_row(nAtomsInRow),stat=alloc_stat)
355      if (alloc_stat /= 0) then
356         status = -1
357         return
# Line 374 | Line 378 | contains
378         status = -1
379         return
380      end if
377 –    if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
378 –    allocate(d2frhodrhodrho_col(nAtomsInCol),stat=alloc_stat)
379 –    if (alloc_stat /= 0) then
380 –       status = -1
381 –       return
382 –    end if
381  
382   #endif
383  
# Line 397 | Line 395 | contains
395  
396    end subroutine setCutoffSC
397  
398 +  subroutine useGeometricMixing()
399 +    useGeometricDistanceMixing = .true.
400 +    haveMixingMap = .false.
401 +    return
402 +  end subroutine useGeometricMixing
403 +  
404  
405  
402 –  subroutine clean_EAM()
406  
407 +
408 +
409 +
410 +
411 +
412 +  subroutine clean_SC()
413 +
414      ! clean non-IS_MPI first
415      frho = 0.0_dp
416      rho  = 0.0_dp
# Line 415 | Line 425 | contains
425      dfrhodrho_row = 0.0_dp
426      dfrhodrho_col = 0.0_dp
427   #endif
428 <  end subroutine clean_EAM
428 >  end subroutine clean_SC
429  
430  
431  
# Line 431 | Line 441 | contains
441      real(kind = dp) :: rho_j_at_i
442  
443      ! we don't use the derivatives, dummy variables
444 <    real( kind = dp) :: drho,d2rho
444 >    real( kind = dp) :: drho
445      integer :: sc_err
446      
447      integer :: atid1,atid2 ! Global atid    
# Line 452 | Line 462 | contains
462      Atid2 = Atid(Atom2)
463   #endif
464  
465 <    Myid_atom1 = SCTypeList%atidtoSCtype(Atid1)
466 <    Myid_atom2 = SCTypeList%atidtoSCtype(Atid2)
465 >    Myid_atom1 = SCList%atidtoSCtype(Atid1)
466 >    Myid_atom2 = SCList%atidtoSCtype(Atid2)
467  
468  
469  
# Line 473 | Line 483 | contains
483  
484  
485   !! Calculate the rho_a for all local atoms
486 <  subroutine calc_eam_preforce_Frho(nlocal,pot)
486 >  subroutine calc_sc_preforce_Frho(nlocal,pot)
487      integer :: nlocal
488      real(kind=dp) :: pot
489      integer :: i,j
490      integer :: atom
491      real(kind=dp) :: U,U1,U2
492      integer :: atype1
493 <    integer :: me,atid1
494 <    integer :: n_rho_points
493 >    integer :: atid1
494 >    integer :: myid
495  
496  
497      cleanme = .true.
# Line 506 | Line 516 | contains
516  
517      !! Calculate F(rho) and derivative
518      do atom = 1, nlocal
519 <       Myid = ScTypeList%atidtoSctype(Atid(atom))
520 <       frho(atom) = -MixingMap(Myid,Myid)%c * MixingMap(Myid,Myid)%epsilon &
521 <            * sqrt(rho(i))
519 >       Myid = SCList%atidtoSctype(Atid(atom))
520 >       frho(atom) = -SCList%SCTypes(Myid)%c * &
521 >            SCList%SCTypes(Myid)%epsilon * sqrt(rho(i))
522 >
523         dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
513 –       d2frhodrhodrho(atom) = -0.5_dp*dfrhodrho(atom)/rho(atom)
524         pot = pot + u
525      enddo
526  
# Line 532 | Line 542 | contains
542      if (sc_err /=  0) then
543         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
544      endif
535 –
536 –    if (nmflag) then
537 –       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_atom_row)
538 –       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_atom_col)
539 –    endif
545   #endif
546      
547      
548 <  end subroutine calc_eam_preforce_Frho
548 >  end subroutine calc_sc_preforce_Frho
549  
550  
551    !! Does Sutton-Chen  pairwise Force calculation.  
552 <  subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
552 >  subroutine do_sc_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
553         pot, f, do_pot)
554      !Arguments    
555      integer, intent(in) ::  atom1, atom2
# Line 557 | Line 562 | contains
562      logical, intent(in) :: do_pot
563  
564      real( kind = dp ) :: drdx,drdy,drdz
565 <    real( kind = dp ) :: d2
566 <    real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
562 <    real( kind = dp ) :: rha,drha,d2rha, dpha
563 <    real( kind = dp ) :: rhb,drhb,d2rhb, dphb
565 >    real( kind = dp ) :: dvpdr
566 >    real( kind = dp ) :: drhodr
567      real( kind = dp ) :: dudr
568      real( kind = dp ) :: rcij
569      real( kind = dp ) :: drhoidr,drhojdr
567 –    real( kind = dp ) :: d2rhoidrdr
568 –    real( kind = dp ) :: d2rhojdrdr
570      real( kind = dp ) :: Fx,Fy,Fz
571      real( kind = dp ) :: r,d2pha,phb,d2phb
572 <
572 >    real( kind = dp ) :: pot_temp,vptmp
573 >    real( kind = dp ) :: epsilonij,aij,nij,mij,vcij
574      integer :: id1,id2
575      integer  :: mytype_atom1
576      integer  :: mytype_atom2
# Line 577 | Line 579 | contains
579  
580      ! write(*,*) "Frho: ", Frho(atom1)
581  
582 <    phab = 0.0E0_DP
582 >
583      dvpdr = 0.0E0_DP
582 –    d2vpdrdr = 0.0E0_DP
584  
585  
586   #ifdef IS_MPI
# Line 590 | Line 591 | contains
591         atid2 = atid(atom2)
592   #endif
593  
594 <       mytype_atom1 = SCTypeList%atidToSCType(atid1)
595 <       mytype_atom2 = SCTypeList%atidTOSCType(atid2)
594 >       mytype_atom1 = SCList%atidToSCType(atid1)
595 >       mytype_atom2 = SCList%atidTOSCType(atid2)
596  
597         drdx = d(1)/rij
598         drdy = d(2)/rij
# Line 604 | Line 605 | contains
605         mij       = MixingMap(mytype_atom1,mytype_atom2)%m
606         vcij      = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot              
607  
608 <       vptmp = dij*((aij/r)**nij)
608 >       vptmp = epsilonij*((aij/r)**nij)
609  
610  
611         dvpdr = -nij*vptmp/r
611 –       d2vpdrdr = -dvpdr*(nij+1.0d0)/r
612 –      
612         drhodr = -mij*((aij/r)**mij)/r
613 <       d2rhodrdr = -drhodr*(mij+1.0d0)/r
613 >
614        
615 <       dudr = drhodr*(dfrhodrho(i)+dfrhodrho(j)) &
615 >       dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
616              + dvpdr
617 +      
618 +       pot_temp = vptmp + vcij
619  
619 –
620  
621   #ifdef IS_MPI
622         dudr = drhodr*(dfrhodrho_row(atom1)+dfrhodrho_col(atom2)) &
# Line 635 | Line 635 | contains
635  
636   #ifdef IS_MPI
637         if (do_pot) then
638 <          pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + phab*0.5
639 <          pot_Col(METALLIC_POT,atom2) = pot_Col(METALLIC_POT,atom2) + phab*0.5
638 >          pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + (pot_temp)*0.5
639 >          pot_Col(METALLIC_POT,atom2) = pot_Col(METALLIC_POT,atom2) + (pot_temp)*0.5
640         end if
641  
642         f_Row(1,atom1) = f_Row(1,atom1) + fx
# Line 649 | Line 649 | contains
649   #else
650  
651         if(do_pot) then
652 <          pot = pot + phab
652 >          pot = pot + pot_temp
653         end if
654  
655         f(1,atom1) = f(1,atom1) + fx
# Line 661 | Line 661 | contains
661         f(3,atom2) = f(3,atom2) - fz
662   #endif
663  
664 <       vpair = vpair + phab
664 >      
665   #ifdef IS_MPI
666         id1 = AtomRowToGlobal(atom1)
667         id2 = AtomColToGlobal(atom2)
# Line 678 | Line 678 | contains
678  
679         endif
680  
681 –       if (nmflag) then
681  
682 <          drhoidr = drha
684 <          drhojdr = drhb
685 <          d2rhoidrdr = d2rha
686 <          d2rhojdrdr = d2rhb
682 >  end subroutine do_sc_pair
683  
688 – #ifdef IS_MPI
689 –          d2 = d2vpdrdr + &
690 –               d2rhoidrdr*dfrhodrho_col(atom2) + &
691 –               d2rhojdrdr*dfrhodrho_row(atom1) + &
692 –               drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
693 –               drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
684  
695 – #else
685  
697 –          d2 = d2vpdrdr + &
698 –               d2rhoidrdr*dfrhodrho(atom2) + &
699 –               d2rhojdrdr*dfrhodrho(atom1) + &
700 –               drhoidr*drhoidr*d2frhodrhodrho(atom2) + &
701 –               drhojdr*drhojdr*d2frhodrhodrho(atom1)
702 – #endif
703 –       end if
704 –
705 –    endif
706 –  end subroutine do_eam_pair
707 –
708 –
709 –
686   end module suttonchen

Diff Legend

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