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 2413 by chuckv, Thu Nov 3 23:06:00 2005 UTC vs.
Revision 2530 by chuckv, Fri Dec 30 00:18:28 2005 UTC

# 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 63 | Line 64 | module suttonchen
64    integer, save :: SC_Mixing_Policy
65    real(kind = dp), save :: SC_rcut
66    logical, save :: haveRcut = .false.
67 <  logical, save,:: haveMixingMap = .false.
67 >  logical, save :: haveMixingMap = .false.
68 >  logical, save :: useGeometricDistanceMixing = .false.
69  
70 +
71 +
72 +
73    character(len = statusMsgSize) :: errMesg
74    integer :: sc_err
75  
# Line 82 | Line 87 | module suttonchen
87       real(kind=dp)     :: n
88       real(kind=dp)     :: alpha
89       real(kind=dp)     :: epsilon
90 +     real(kind=dp)     :: sc_rcut
91    end type SCtype
92  
93  
# Line 89 | Line 95 | module suttonchen
95    real( kind = dp), dimension(:), allocatable :: rho
96    real( kind = dp), dimension(:), allocatable :: frho
97    real( kind = dp), dimension(:), allocatable :: dfrhodrho
92 –  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
98  
99  
100 +
101    !! Arrays for MPI storage
102   #ifdef IS_MPI
103    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
# Line 101 | Line 107 | module suttonchen
107    real( kind = dp),save, dimension(:), allocatable :: rho_row
108    real( kind = dp),save, dimension(:), allocatable :: rho_col
109    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
110   #endif
111  
112    type, private :: SCTypeList
# Line 121 | Line 125 | module suttonchen
125   type :: MixParameters
126       real(kind=DP) :: alpha
127       real(kind=DP) :: epsilon
128 <     real(kind=dp) :: sigma6
128 >     real(kind=DP) :: m
129 >     real(Kind=DP) :: n
130 >     real(kind=DP) :: vpair_pot
131       real(kind=dp) :: rCut
126 –     real(kind=dp) :: vpair_pot
132       logical       :: rCutWasSet = .false.
133 <     logical       :: shiftedPot
129 <     logical       :: isSoftCore = .false.
133 >
134    end type MixParameters
135  
136    type(MixParameters), dimension(:,:), allocatable :: MixingMap
137  
138  
139  
136 –  public :: init_SC_FF
140    public :: setCutoffSC
141    public :: do_SC_pair
142    public :: newSCtype
143    public :: calc_SC_prepair_rho
144 +  public :: calc_SC_preforce_Frho
145    public :: clean_SC
146    public :: destroySCtypes
147    public :: getSCCut
148 <  public :: setSCDefaultCutoff
149 <  public :: setSCUniformCutoff
148 > ! public :: setSCDefaultCutoff
149 > ! public :: setSCUniformCutoff
150 >
151  
152   contains
153  
154  
155 <  subroutine newSCtype(c,m,n,alpha,epsilon,c_ident,status)
156 <    real (kind = dp )                      :: lattice_constant
157 <    integer                                :: eam_nrho
158 <    real (kind = dp )                      :: eam_drho
159 <    integer                                :: eam_nr
160 <    real (kind = dp )                      :: eam_dr
161 <    real (kind = dp )                      :: rcut
162 <    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
155 >  subroutine newSCtype(c_ident,c,m,n,alpha,epsilon,status)
156 >    real (kind = dp )                      :: c ! Density Scaling
157 >    real (kind = dp )                      :: m ! Density Exponent
158 >    real (kind = dp )                      :: n ! Pair Potential Exponent
159 >    real (kind = dp )                      :: alpha ! Length Scaling
160 >    real (kind = dp )                      :: epsilon ! Energy Scaling
161 >
162 >
163      integer                                :: c_ident
164      integer                                :: status
165  
# Line 170 | Line 173 | contains
173  
174  
175      !! 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.
176  
177  
178      ! check to see if this is the first time into
179 <    if (.not.associated(SCTypeList%EAMParams)) then
180 <       call getMatchingElementList(atypes, "is_SuttonChen", .true., nSCtypes, MatchList)
181 <       SCTypeList%nSCtypes = nSCtypes
182 <       allocate(SCTypeList%SCTypes(nSCTypes))
179 >    if (.not.associated(SCList%SCTypes)) then
180 >       call getMatchingElementList(atypes, "is_SC", .true., nSCtypes, MatchList)
181 >       SCList%nSCtypes = nSCtypes
182 >       allocate(SCList%SCTypes(nSCTypes))
183         nAtypes = getSize(atypes)
184 <       allocate(SCTypeList%atidToSCType(nAtypes))
184 >       allocate(SCList%atidToSCType(nAtypes))
185      end if
186  
187 <    SCTypeList%currentSCType = SCTypeList%currentSCType + 1
188 <    current = SCTypeList%currentSCType
187 >    SCList%currentSCType = SCList%currentSCType + 1
188 >    current = SCList%currentSCType
189  
190      myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
191 <    SCTypeList%atidToSCType(myATID) = current
191 >    SCList%atidToSCType(myATID) = current
192  
193  
194    
195 <    SCTypeList%SCTypes(current)%atid         = c_ident
196 <    SCTypeList%SCTypes(current)%alpha        = alpha
197 <    SCTypeList%SCTypes(current)%c            = c
198 <    SCTypeList%SCTypes(current)%m            = m
199 <    SCTypeList%SCTypes(current)%n            = n
200 <    SCTypeList%SCTypes(current)%epsilon      = epsilon
195 >    SCList%SCTypes(current)%atid         = c_ident
196 >    SCList%SCTypes(current)%alpha        = alpha
197 >    SCList%SCTypes(current)%c            = c
198 >    SCList%SCTypes(current)%m            = m
199 >    SCList%SCTypes(current)%n            = n
200 >    SCList%SCTypes(current)%epsilon      = epsilon
201    end subroutine newSCtype
202  
203    
204    subroutine destroySCTypes()
205 +    if (associated(SCList%SCtypes)) then
206 +       deallocate(SCList%SCTypes)
207 +       SCList%SCTypes=>null()
208 +    end if
209 +    if (associated(SCList%atidToSCtype)) then
210 +       deallocate(SCList%atidToSCtype)
211 +       SCList%atidToSCtype=>null()
212 +    end if
213  
204 –    
205 –    if(allocated(SCTypeList)) deallocate(SCTypeList)
214  
207 –
208 –
215    end subroutine destroySCTypes
216  
217  
218  
219    function getSCCut(atomID) result(cutValue)
220      integer, intent(in) :: atomID
221 <    integer :: eamID
221 >    integer :: scID
222      real(kind=dp) :: cutValue
223      
224 <    eamID = SCTypeList%atidToEAMType(atomID)
225 <    cutValue = SCTypeList%EAMParams(eamID)%eam_rcut
224 >    scID = SCList%atidToSCType(atomID)
225 >    cutValue = 2.0_dp * SCList%SCTypes(scID)%alpha
226    end function getSCCut
227  
228  
# Line 238 | Line 244 | contains
244      if (.not. allocated(MixingMap)) then
245         allocate(MixingMap(nSCtypes, nSCtypes))
246      endif
247 <
247 >    useGeometricDistanceMixing = usesGeometricDistanceMixing()
248      do i = 1, nSCtypes
249  
250         e1 = SCList%SCtypes(i)%epsilon
# Line 254 | Line 260 | contains
260            n2 = SCList%SCtypes(j)%n
261            alpha2 = SCList%SCtypes(j)%alpha
262  
263 <          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
263 >          if (useGeometricDistanceMixing) then
264 >             MixingMap(i,j)%alpha = sqrt(alpha1 * alpha2) !SC formulation
265 >          else
266 >             MixingMap(i,j)%alpha = 0.5_dp * (alpha1 + alpha2) ! Goddard formulation
267 >          endif
268 >
269 >          MixingMap(i,j)%epsilon = sqrt(e1 * e2)
270            MixingMap(i,j)%m = 0.5_dp*(m1+m2)
271            MixingMap(i,j)%n = 0.5_dp*(n1+n2)
272            MixingMap(i,j)%alpha = 0.5_dp*(alpha1+alpha2)
273            MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
274 <          MixingMap(i,j)%vpair_rcut = MixingMap%epsilon(i,j)* &
275 <               (MixingMap%alpha(i,j)/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
276 <
274 >          MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* &
275 >               (MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
276 >          if (i.ne.j) then
277 >             MixingMap(j,i)%epsilon    = MixingMap(i,j)%epsilon
278 >             MixingMap(j,i)%m          = MixingMap(i,j)%m
279 >             MixingMap(j,i)%n          = MixingMap(i,j)%n
280 >             MixingMap(j,i)%alpha      = MixingMap(i,j)%alpha
281 >             MixingMap(j,i)%rcut = MixingMap(i,j)%rcut
282 >             MixingMap(j,i)%vpair_pot = MixingMap(i,j)%vpair_pot
283 >          endif
284         enddo
285      enddo
286      
# Line 311 | Line 330 | contains
330         return
331      end if
332  
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 –
333   #ifdef IS_MPI
334  
335      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 342 | Line 354 | contains
354      end if
355      if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
356      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)
357      if (alloc_stat /= 0) then
358         status = -1
359         return
# Line 374 | Line 380 | contains
380         status = -1
381         return
382      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
383  
384   #endif
385  
# Line 397 | Line 397 | contains
397  
398    end subroutine setCutoffSC
399  
400 +  subroutine clean_SC()
401  
401 –
402 –  subroutine clean_EAM()
403 –
402      ! clean non-IS_MPI first
403      frho = 0.0_dp
404      rho  = 0.0_dp
# Line 415 | Line 413 | contains
413      dfrhodrho_row = 0.0_dp
414      dfrhodrho_col = 0.0_dp
415   #endif
416 <  end subroutine clean_EAM
416 >  end subroutine clean_SC
417  
418  
419  
420    !! Calculates rho_r
421 <  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq)
421 >  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq, rcut)
422      integer :: atom1,atom2
423      real(kind = dp), dimension(3) :: d
424      real(kind = dp), intent(inout)               :: r
425 <    real(kind = dp), intent(inout)               :: rijsq
425 >    real(kind = dp), intent(inout)               :: rijsq, rcut
426      ! value of electron density rho do to atom i at atom j
427      real(kind = dp) :: rho_i_at_j
428      ! value of electron density rho do to atom j at atom i
429      real(kind = dp) :: rho_j_at_i
430  
431      ! we don't use the derivatives, dummy variables
432 <    real( kind = dp) :: drho,d2rho
432 >    real( kind = dp) :: drho
433      integer :: sc_err
434      
435      integer :: atid1,atid2 ! Global atid    
# Line 452 | Line 450 | contains
450      Atid2 = Atid(Atom2)
451   #endif
452  
453 <    Myid_atom1 = SCTypeList%atidtoSCtype(Atid1)
454 <    Myid_atom2 = SCTypeList%atidtoSCtype(Atid2)
453 >    Myid_atom1 = SCList%atidtoSCtype(Atid1)
454 >    Myid_atom2 = SCList%atidtoSCtype(Atid2)
455  
456  
457  
# Line 473 | Line 471 | contains
471  
472  
473   !! Calculate the rho_a for all local atoms
474 <  subroutine calc_eam_preforce_Frho(nlocal,pot)
474 >  subroutine calc_sc_preforce_Frho(nlocal,pot)
475      integer :: nlocal
476      real(kind=dp) :: pot
477      integer :: i,j
478      integer :: atom
479      real(kind=dp) :: U,U1,U2
480      integer :: atype1
481 <    integer :: me,atid1
482 <    integer :: n_rho_points
481 >    integer :: atid1
482 >    integer :: myid
483  
484  
485      cleanme = .true.
# Line 506 | Line 504 | contains
504  
505      !! Calculate F(rho) and derivative
506      do atom = 1, nlocal
507 <       Myid = ScTypeList%atidtoSctype(Atid(atom))
508 <       frho(atom) = -MixingMap(Myid,Myid)%c * MixingMap(Myid,Myid)%epsilon &
509 <            * sqrt(rho(i))
507 >       Myid = SCList%atidtoSctype(Atid(atom))
508 >       frho(atom) = -SCList%SCTypes(Myid)%c * &
509 >            SCList%SCTypes(Myid)%epsilon * sqrt(rho(i))
510 >
511         dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
513 –       d2frhodrhodrho(atom) = -0.5_dp*dfrhodrho(atom)/rho(atom)
512         pot = pot + u
513      enddo
514  
515 <    #ifdef IS_MPI
515 > #ifdef IS_MPI
516      !! communicate f(rho) and derivatives back into row and column arrays
517      call gather(frho,frho_row,plan_atom_row, sc_err)
518      if (sc_err /=  0) then
# Line 532 | Line 530 | contains
530      if (sc_err /=  0) then
531         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
532      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
533   #endif
534      
535      
536 <  end subroutine calc_eam_preforce_Frho
536 >  end subroutine calc_sc_preforce_Frho
537  
538  
539    !! Does Sutton-Chen  pairwise Force calculation.  
540 <  subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
540 >  subroutine do_sc_pair(atom1, atom2, d, rij, r2, rcut, sw, vpair, fpair, &
541         pot, f, do_pot)
542      !Arguments    
543      integer, intent(in) ::  atom1, atom2
544 <    real( kind = dp ), intent(in) :: rij, r2
544 >    real( kind = dp ), intent(in) :: rij, r2, rcut
545      real( kind = dp ) :: pot, sw, vpair
546      real( kind = dp ), dimension(3,nLocal) :: f
547      real( kind = dp ), intent(in), dimension(3) :: d
# Line 557 | Line 550 | contains
550      logical, intent(in) :: do_pot
551  
552      real( kind = dp ) :: drdx,drdy,drdz
553 <    real( kind = dp ) :: d2
554 <    real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
562 <    real( kind = dp ) :: rha,drha,d2rha, dpha
563 <    real( kind = dp ) :: rhb,drhb,d2rhb, dphb
553 >    real( kind = dp ) :: dvpdr
554 >    real( kind = dp ) :: drhodr
555      real( kind = dp ) :: dudr
556      real( kind = dp ) :: rcij
557      real( kind = dp ) :: drhoidr,drhojdr
567 –    real( kind = dp ) :: d2rhoidrdr
568 –    real( kind = dp ) :: d2rhojdrdr
558      real( kind = dp ) :: Fx,Fy,Fz
559      real( kind = dp ) :: r,d2pha,phb,d2phb
560 <
560 >    real( kind = dp ) :: pot_temp,vptmp
561 >    real( kind = dp ) :: epsilonij,aij,nij,mij,vcij
562      integer :: id1,id2
563      integer  :: mytype_atom1
564      integer  :: mytype_atom2
# Line 577 | Line 567 | contains
567  
568      ! write(*,*) "Frho: ", Frho(atom1)
569  
570 <    phab = 0.0E0_DP
570 >
571      dvpdr = 0.0E0_DP
582 –    d2vpdrdr = 0.0E0_DP
572  
573  
574   #ifdef IS_MPI
# Line 590 | Line 579 | contains
579         atid2 = atid(atom2)
580   #endif
581  
582 <       mytype_atom1 = SCTypeList%atidToSCType(atid1)
583 <       mytype_atom2 = SCTypeList%atidTOSCType(atid2)
582 >       mytype_atom1 = SCList%atidToSCType(atid1)
583 >       mytype_atom2 = SCList%atidTOSCType(atid2)
584  
585         drdx = d(1)/rij
586         drdy = d(2)/rij
# Line 604 | Line 593 | contains
593         mij       = MixingMap(mytype_atom1,mytype_atom2)%m
594         vcij      = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot              
595  
596 <       vptmp = dij*((aij/r)**nij)
596 >       vptmp = epsilonij*((aij/r)**nij)
597  
598  
599         dvpdr = -nij*vptmp/r
611 –       d2vpdrdr = -dvpdr*(nij+1.0d0)/r
612 –      
600         drhodr = -mij*((aij/r)**mij)/r
601 <       d2rhodrdr = -drhodr*(mij+1.0d0)/r
601 >
602        
603 <       dudr = drhodr*(dfrhodrho(i)+dfrhodrho(j)) &
603 >       dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
604              + dvpdr
605 +      
606 +       pot_temp = vptmp + vcij
607  
619 –
608  
609   #ifdef IS_MPI
610         dudr = drhodr*(dfrhodrho_row(atom1)+dfrhodrho_col(atom2)) &
# Line 635 | Line 623 | contains
623  
624   #ifdef IS_MPI
625         if (do_pot) then
626 <          pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + phab*0.5
627 <          pot_Col(METALLIC_POT,atom2) = pot_Col(METALLIC_POT,atom2) + phab*0.5
626 >          pot_Row(METALLIC_POT,atom1) = pot_Row(METALLIC_POT,atom1) + (pot_temp)*0.5
627 >          pot_Col(METALLIC_POT,atom2) = pot_Col(METALLIC_POT,atom2) + (pot_temp)*0.5
628         end if
629  
630         f_Row(1,atom1) = f_Row(1,atom1) + fx
# Line 649 | Line 637 | contains
637   #else
638  
639         if(do_pot) then
640 <          pot = pot + phab
640 >          pot = pot + pot_temp
641         end if
642  
643         f(1,atom1) = f(1,atom1) + fx
# Line 661 | Line 649 | contains
649         f(3,atom2) = f(3,atom2) - fz
650   #endif
651  
652 <       vpair = vpair + phab
652 >      
653   #ifdef IS_MPI
654         id1 = AtomRowToGlobal(atom1)
655         id2 = AtomColToGlobal(atom2)
# Line 678 | Line 666 | contains
666  
667         endif
668  
681 –       if (nmflag) then
669  
670 <          drhoidr = drha
684 <          drhojdr = drhb
685 <          d2rhoidrdr = d2rha
686 <          d2rhojdrdr = d2rhb
670 >  end subroutine do_sc_pair
671  
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)
672  
695 – #else
673  
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 –
674   end module suttonchen

Diff Legend

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