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 2406 by chuckv, Tue Nov 1 23:32:25 2005 UTC vs.
Revision 2412 by chuckv, Wed Nov 2 23:50:56 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.
67  
68    character(len = statusMsgSize) :: errMesg
69    integer :: eam_err
# Line 85 | Line 86 | module suttonchen
86  
87  
88    !! Arrays for derivatives used in force calculation
88 –  real( kind = dp), dimension(:), allocatable :: frho
89    real( kind = dp), dimension(:), allocatable :: rho
90 <
90 >  real( kind = dp), dimension(:), allocatable :: frho
91    real( kind = dp), dimension(:), allocatable :: dfrhodrho
92    real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
93  
# Line 123 | Line 123 | module suttonchen
123       real(kind=DP) :: epsilon
124       real(kind=dp) :: sigma6
125       real(kind=dp) :: rCut
126 <     real(kind=dp) :: delta
126 >     real(kind=dp) :: vpair_pot
127       logical       :: rCutWasSet = .false.
128       logical       :: shiftedPot
129       logical       :: isSoftCore = .false.
# Line 141 | Line 141 | module suttonchen
141    public :: clean_SC
142    public :: destroySCtypes
143    public :: getSCCut
144 +  public :: setSCDefaultCutoff
145 +  public :: setSCUniformCutoff
146  
147   contains
148  
# Line 158 | Line 160 | contains
160      integer                                :: c_ident
161      integer                                :: status
162  
163 <    integer                                :: nAtypes,nEAMTypes,myATID
163 >    integer                                :: nAtypes,nSCTypes,myATID
164      integer                                :: maxVals
165      integer                                :: alloc_stat
166      integer                                :: current
# Line 173 | Line 175 | contains
175  
176      ! check to see if this is the first time into
177      if (.not.associated(EAMList%EAMParams)) then
178 <       call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
179 <       EAMList%n_eam_types = nEAMtypes
180 <       allocate(EAMList%EAMParams(nEAMTypes))
178 >       call getMatchingElementList(atypes, "is_SuttonChen", .true., nSCtypes, MatchList)
179 >       SCTypeList%nSCtypes = nSCtypes
180 >       allocate(SCTypeList%SCTypes(nSCTypes))
181         nAtypes = getSize(atypes)
182 <       allocate(EAMList%atidToEAMType(nAtypes))
182 >       allocate(SCTypeList%atidToSCType(nAtypes))
183      end if
184  
185 <    EAMList%currentAddition = EAMList%currentAddition + 1
186 <    current = EAMList%currentAddition
185 >    SCTypeList%currentSCType = SCTypeList%currentSCType + 1
186 >    current = SCTypeList%currentSCType
187  
188      myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
189 <    EAMList%atidToEAMType(myATID) = current
189 >    SCTypeList%atidToSCType(myATID) = current
190  
189 –    call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
190 –    if (alloc_stat /= 0) then
191 –       status = -1
192 –       return
193 –    end if
191  
192    
193 <    EAMList%EAMParams(current)%eam_atype    = c_ident
194 <    EAMList%EAMParams(current)%eam_lattice  = lattice_constant
195 <    EAMList%EAMParams(current)%eam_nrho     = eam_nrho
196 <    EAMList%EAMParams(current)%eam_drho     = eam_drho
197 <    EAMList%EAMParams(current)%eam_nr       = eam_nr
198 <    EAMList%EAMParams(current)%eam_dr       = eam_dr
199 <    EAMList%EAMParams(current)%eam_rcut     = rcut
203 <    EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
204 <    EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
205 <    EAMList%EAMParams(current)%eam_F_rho    = eam_F_rho
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
199 >  end subroutine newSCtype
200  
201 <  end subroutine newEAMtype
201 >  
202 >  subroutine destroySCTypes()
203  
204 +    
205 +    if(allocated(SCTypeList)) deallocate(SCTypeList)
206  
210 –  ! kills all eam types entered and sets simulation to uninitalized
211 –  subroutine destroyEAMtypes()
212 –    integer :: i
213 –    type(EAMType), pointer :: tempEAMType=>null()
207  
215 –    do i = 1, EAMList%n_eam_types
216 –       tempEAMType => eamList%EAMParams(i)
217 –       call deallocate_EAMType(tempEAMType)
218 –    end do
219 –    if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
220 –    eamList%EAMParams => null()
208  
209 <    eamList%n_eam_types = 0
223 <    eamList%currentAddition = 0
209 >  end subroutine destroySCTypes
210  
225 –  end subroutine destroyEAMtypes
211  
212 <  function getEAMCut(atomID) result(cutValue)
212 >
213 >  function getSCCut(atomID) result(cutValue)
214      integer, intent(in) :: atomID
215      integer :: eamID
216      real(kind=dp) :: cutValue
217      
218      eamID = EAMList%atidToEAMType(atomID)
219      cutValue = EAMList%EAMParams(eamID)%eam_rcut
220 <  end function getEAMCut
220 >  end function getSCCut
221  
236 –  subroutine init_EAM_FF(status)
237 –    integer :: status
238 –    integer :: i,j
239 –    real(kind=dp) :: current_rcut_max
240 –    integer :: alloc_stat
241 –    integer :: number_r, number_rho
222  
223  
224 <    status = 0
225 <    if (EAMList%currentAddition == 0) then
226 <       call handleError("init_EAM_FF","No members in EAMList")
227 <       status = -1
224 >
225 >  subroutine createMixingMap()
226 >    integer :: nSCtypes, i, j
227 >    real ( kind = dp ) :: e1, e2,m1,m2,alpha1,alpha2,n1,n2
228 >    real ( kind = dp ) :: rcut6, tp6, tp12
229 >    logical :: isSoftCore1, isSoftCore2, doShift
230 >
231 >    if (SCList%currentSCtype == 0) then
232 >       call handleError("SuttonChen", "No members in SCMap")
233         return
234      end if
235  
236 +    nSCtypes = SCList%nSCtypes
237  
238 <    do i = 1, EAMList%currentAddition
238 >    if (.not. allocated(MixingMap)) then
239 >       allocate(MixingMap(nSCtypes, nSCtypes))
240 >    endif
241  
242 <       ! Build array of r values
242 >    do i = 1, nSCtypes
243  
244 <       do j = 1,EAMList%EAMParams(i)%eam_nr
245 <          EAMList%EAMParams(i)%eam_rvals(j) = &
246 <               real(j-1,kind=dp)* &
247 <               EAMList%EAMParams(i)%eam_dr
260 <       end do
261 <       ! Build array of rho values
262 <       do j = 1,EAMList%EAMParams(i)%eam_nrho
263 <          EAMList%EAMParams(i)%eam_rhovals(j) = &
264 <               real(j-1,kind=dp)* &
265 <               EAMList%EAMParams(i)%eam_drho
266 <       end do
267 <       ! convert from eV to kcal / mol:
268 <       EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
244 >       e1 = SCList%SCtypes(i)%epsilon
245 >       m1 = SCList%SCtypes(i)%m
246 >       n1 = SCList%SCtypes(i)%n
247 >       alpha1 = SCList%SCtypes(i)%alpha
248  
249 <       ! precompute the pair potential and get it into kcal / mol:
250 <       EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
272 <       do j = 2, EAMList%EAMParams(i)%eam_nr
273 <          EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
274 <          EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
275 <       enddo
276 <    end do
249 >       do j = i, nSCtypes
250 >          
251  
252 <
253 <    do i = 1,  EAMList%currentAddition
254 <       number_r   = EAMList%EAMParams(i)%eam_nr
255 <       number_rho = EAMList%EAMParams(i)%eam_nrho
252 >          e2 = SCList%SCtypes(j)%epsilon
253 >          m2 = SCList%SCtypes(j)%m
254 >          n2 = SCList%SCtypes(j)%n
255 >          alpha2 = SCList%SCtypes(j)%alpha
256  
257 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
258 <            EAMList%EAMParams(i)%eam_rho_r, &
259 <            EAMList%EAMParams(i)%eam_rho_r_pp, &
260 <            0.0E0_DP, 0.0E0_DP, 'N')
261 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
262 <            EAMList%EAMParams(i)%eam_Z_r, &
263 <            EAMList%EAMParams(i)%eam_Z_r_pp, &
264 <            0.0E0_DP, 0.0E0_DP, 'N')
265 <       call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
292 <            EAMList%EAMParams(i)%eam_F_rho, &
293 <            EAMList%EAMParams(i)%eam_F_rho_pp, &
294 <            0.0E0_DP, 0.0E0_DP, 'N')
295 <       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
296 <            EAMList%EAMParams(i)%eam_phi_r, &
297 <            EAMList%EAMParams(i)%eam_phi_r_pp, &
298 <            0.0E0_DP, 0.0E0_DP, 'N')
257 >          MixingMap(i,j)%epsilon = dsqrt(e1 * e2)
258 >          MixingMap(i,j)%m = 0.5_dp*(m1+m2)
259 >          MixingMap(i,j)%n = 0.5_dp*(n1+n2)
260 >          MixingMap(i,j)%alpha = 0.5_dp*(alpha1+alpha2)
261 >          MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha
262 >          MixingMap(i,j)%vpair_rcut = MixingMap%epsilon(i,j)* &
263 >               (MixingMap%alpha(i,j)/MixingMap(i,j)%rcut)**MixingMap(i,j)%n
264 >
265 >       enddo
266      enddo
267 +    
268 +    haveMixingMap = .true.
269 +    
270 +  end subroutine createMixingMap
271 +  
272  
301 –    !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
302 –    !! find the smallest rcut for any eam atype
303 –    !       do i = 2, EAMList%currentAddition
304 –    !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
305 –    !       end do
306 –
307 –    !       EAM_rcut = current_rcut_max
308 –    !       EAM_rcut_orig = current_rcut_max
309 –    !       do i = 1, EAMList%currentAddition
310 –    !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
311 –    !       end do
312 –    !! Allocate arrays for force calculation
313 –
314 –    call allocateEAM(alloc_stat)
315 –    if (alloc_stat /= 0 ) then
316 –       write(*,*) "allocateEAM failed"
317 –       status = -1
318 –       return
319 –    endif
320 –
321 –  end subroutine init_EAM_FF
322 –
273    !! routine checks to see if array is allocated, deallocates array if allocated
274    !! and then creates the array to the required size
275 <  subroutine allocateEAM(status)
275 >  subroutine allocateSC(status)
276      integer, intent(out) :: status
277  
278   #ifdef IS_MPI
# Line 338 | Line 288 | contains
288      nAtomsInCol = getNatomsInCol(plan_atom_col)
289   #endif
290  
291 +
292 +
293      if (allocated(frho)) deallocate(frho)
294      allocate(frho(nlocal),stat=alloc_stat)
295 <    if (alloc_stat /= 0) then
295 >    if (alloc_stat /= 0) then
296         status = -1
297         return
298      end if
299 +
300      if (allocated(rho)) deallocate(rho)
301      allocate(rho(nlocal),stat=alloc_stat)
302      if (alloc_stat /= 0) then
# Line 430 | Line 383 | contains
383  
384   #endif
385  
386 <  end subroutine allocateEAM
386 >  end subroutine allocateSC
387  
388    !! C sets rcut to be the largest cutoff of any atype
389    !! present in this simulation. Doesn't include all atypes
390    !! sim knows about, just those in the simulation.
391 <  subroutine setCutoffEAM(rcut, status)
391 >  subroutine setCutoffSC(rcut, status)
392      real(kind=dp) :: rcut
393      integer :: status
394      status = 0
395  
396 <    EAM_rcut = rcut
396 >    SC_rcut = rcut
397  
398 <  end subroutine setCutoffEAM
398 >  end subroutine setCutoffSC
399  
400  
401  
# Line 463 | Line 416 | contains
416      dfrhodrho_col = 0.0_dp
417   #endif
418    end subroutine clean_EAM
466 –
467 –
468 –
469 –  subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
470 –    integer, intent(in)          :: eam_n_rho
471 –    integer, intent(in)          :: eam_n_r
472 –    type (EAMType)               :: thisEAMType
473 –    integer, optional   :: stat
474 –    integer             :: alloc_stat
475 –
476 –
477 –
478 –    if (present(stat)) stat = 0
479 –
480 –    allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
481 –    if (alloc_stat /= 0 ) then
482 –       if (present(stat)) stat = -1
483 –       return
484 –    end if
485 –    allocate(thisEAMType%eam_rhovals(eam_n_rho),stat=alloc_stat)  
486 –    if (alloc_stat /= 0 ) then
487 –       if (present(stat)) stat = -1
488 –       return
489 –    end if
490 –    allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat)  
491 –    if (alloc_stat /= 0 ) then
492 –       if (present(stat)) stat = -1
493 –       return
494 –    end if
495 –    allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat)        
496 –    if (alloc_stat /= 0 ) then
497 –       if (present(stat)) stat = -1
498 –       return
499 –    end if
500 –    allocate(thisEAMType%eam_rho_r(eam_n_r),stat=alloc_stat)      
501 –    if (alloc_stat /= 0 ) then
502 –       if (present(stat)) stat = -1
503 –       return
504 –    end if
505 –    allocate(thisEAMType%eam_phi_r(eam_n_r),stat=alloc_stat)      
506 –    if (alloc_stat /= 0 ) then
507 –       if (present(stat)) stat = -1
508 –       return
509 –    end if
510 –    allocate(thisEAMType%eam_F_rho_pp(eam_n_rho),stat=alloc_stat)  
511 –    if (alloc_stat /= 0 ) then
512 –       if (present(stat)) stat = -1
513 –       return
514 –    end if
515 –    allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat)  
516 –    if (alloc_stat /= 0 ) then
517 –       if (present(stat)) stat = -1
518 –       return
519 –    end if
520 –    allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat)  
521 –    if (alloc_stat /= 0 ) then
522 –       if (present(stat)) stat = -1
523 –       return
524 –    end if
525 –    allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat)
526 –    if (alloc_stat /= 0 ) then
527 –       if (present(stat)) stat = -1
528 –       return
529 –    end if
419  
420  
532 –  end subroutine allocate_EAMType
421  
534 –
535 –  subroutine deallocate_EAMType(thisEAMType)
536 –    type (EAMtype), pointer :: thisEAMType
537 –
538 –    ! free Arrays in reverse order of allocation...
539 –    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
540 –    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
541 –    if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)  
542 –    if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)  
543 –    if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)      
544 –    if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)      
545 –    if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)  
546 –    if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
547 –    if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
548 –    if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
549 –
550 –  end subroutine deallocate_EAMType
551 –
422    !! Calculates rho_r
423 <  subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
423 >  subroutine calc_sc_prepair_rho(atom1,atom2,d,r,rijsq)
424      integer :: atom1,atom2
425      real(kind = dp), dimension(3) :: d
426      real(kind = dp), intent(inout)               :: r
# Line 582 | Line 452 | contains
452      Atid2 = Atid(Atom2)
453   #endif
454  
455 <    Myid_atom1 = Eamlist%atidtoeamtype(Atid1)
456 <    Myid_atom2 = Eamlist%atidtoeamtype(Atid2)
455 >    Myid_atom1 = SCTypeList%atidtoSCtype(Atid1)
456 >    Myid_atom2 = SCTypeList%atidtoSCtype(Atid2)
457  
588 –    if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
458  
459  
460 +       rho_i_at_j = (MixingMap(Myid_atom1,Myid_atom2)%alpha/r)&
461 +            **MixingMap(Myid_atom1,Myid_atom2)%m
462 +       rho_j_at_i = rho_i_at_j
463  
592 –       call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
593 –            EAMList%EAMParams(myid_atom1)%eam_rvals, &
594 –            EAMList%EAMParams(myid_atom1)%eam_rho_r, &
595 –            EAMList%EAMParams(myid_atom1)%eam_rho_r_pp, &
596 –            r, rho_i_at_j,drho,d2rho)
597 –
598 –
599 –
464   #ifdef  IS_MPI
465         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
602 – #else
603 –       rho(atom2) = rho(atom2) + rho_i_at_j
604 – #endif
605 –             ! write(*,*) atom1,atom2,r,rho_i_at_j
606 –    endif
607 –
608 –    if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
609 –       call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
610 –            EAMList%EAMParams(myid_atom2)%eam_rvals, &
611 –            EAMList%EAMParams(myid_atom2)%eam_rho_r, &
612 –            EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
613 –            r, rho_j_at_i,drho,d2rho)
614 –
615 –
616 –
617 –
618 – #ifdef  IS_MPI
466         rho_row(atom1) = rho_row(atom1) + rho_j_at_i
467   #else
468 +       rho(atom2) = rho(atom2) + rho_i_at_j
469         rho(atom1) = rho(atom1) + rho_j_at_i
470   #endif
623 –    endif
471  
472  
473  
474 +  end subroutine calc_sc_prepair_rho
475  
476  
477 <
630 <  end subroutine calc_eam_prepair_rho
631 <
632 <
633 <
634 <
635 <  !! Calculate the functional F(rho) for all local atoms
477 > !! Calculate the functional F(rho) for all local atoms
478    subroutine calc_eam_preforce_Frho(nlocal,pot)
479      integer :: nlocal
480      real(kind=dp) :: pot
# Line 656 | Line 498 | contains
498      if (eam_err /= 0 ) then
499         write(errMsg,*) " Error scattering rho_col into rho"
500         call handleError(RoutineName,errMesg)
501 +
502      endif
503  
504      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
# Line 663 | Line 506 | contains
506  
507  
508  
509 <    !! Calculate F(rho) and derivative
509 >    !! Calculate F(rho) and derivative
510      do atom = 1, nlocal
511 <       atid1 = atid(atom)
512 <       me = eamList%atidToEAMtype(atid1)
513 <       n_rho_points = EAMList%EAMParams(me)%eam_nrho
514 <       !  Check to see that the density is not greater than the larges rho we have calculated
515 <       if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
673 <          call eam_splint(n_rho_points, &
674 <               EAMList%EAMParams(me)%eam_rhovals, &
675 <               EAMList%EAMParams(me)%eam_f_rho, &
676 <               EAMList%EAMParams(me)%eam_f_rho_pp, &
677 <               rho(atom), & ! Actual Rho
678 <               u, u1, u2)
679 <       else
680 <          ! Calculate F(rho with the largest available rho value
681 <          call eam_splint(n_rho_points, &
682 <               EAMList%EAMParams(me)%eam_rhovals, &
683 <               EAMList%EAMParams(me)%eam_f_rho, &
684 <               EAMList%EAMParams(me)%eam_f_rho_pp, &
685 <               EAMList%EAMParams(me)%eam_rhovals(n_rho_points), & ! Largest rho
686 <               u,u1,u2)
687 <       end if
688 <
689 <
690 <       frho(atom) = u
691 <       dfrhodrho(atom) = u1
692 <       d2frhodrhodrho(atom) = u2
511 >       Myid = ScTypeList%atidtoSctype(Atid(atom))
512 >       frho(atom) = -MixingMap(Myid,Myid)%c * MixingMap(Myid,Myid)%epsilon &
513 >            * sqrt(rho(i))
514 >       dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom)
515 >       d2frhodrhodrho(atom) = -0.5_dp*dfrhodrho(atom)/rho(atom)
516         pot = pot + u
694 –
517      enddo
518  
519 <
698 <
699 < #ifdef IS_MPI
519 >    #ifdef IS_MPI
520      !! communicate f(rho) and derivatives back into row and column arrays
521      call gather(frho,frho_row,plan_atom_row, eam_err)
522      if (eam_err /=  0) then
# Line 715 | Line 535 | contains
535         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
536      endif
537  
718 –
719 –
720 –
721 –
538      if (nmflag) then
539         call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_atom_row)
540         call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_atom_col)
541      endif
542   #endif
543 <
544 <
543 >    
544 >    
545    end subroutine calc_eam_preforce_Frho
546  
547  
548 <
733 <
734 <  !! Does EAM pairwise Force calculation.  
548 >  !! Does Sutton-Chen  pairwise Force calculation.  
549    subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
550         pot, f, do_pot)
551      !Arguments    
# Line 935 | Line 749 | contains
749  
750      endif
751    end subroutine do_eam_pair
938 –
939 –
940 –  subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
941 –
942 –    integer :: atype, nx, j
943 –    real( kind = DP ), dimension(:) :: xa
944 –    real( kind = DP ), dimension(:) :: ya
945 –    real( kind = DP ), dimension(:) :: yppa
946 –    real( kind = DP ) :: x, y
947 –    real( kind = DP ) :: dy, d2y
948 –    real( kind = DP ) :: del, h, a, b, c, d
949 –    integer :: pp_arraySize
950 –
951 –
952 –    ! this spline code assumes that the x points are equally spaced
953 –    ! do not attempt to use this code if they are not.
954 –
955 –
956 –    ! find the closest point with a value below our own:
957 –    j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
958 –
959 –    ! check to make sure we're inside the spline range:
960 –    if ((j.gt.nx).or.(j.lt.1)) then
961 –       write(errMSG,*) "EAM_splint: x is outside bounds of spline: ",x,j
962 –       call handleError(routineName,errMSG)
963 –    endif
964 –    ! check to make sure we haven't screwed up the calculation of j:
965 –    if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
966 –       if (j.ne.nx) then
967 –          write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
968 –          call handleError(routineName,errMSG)
969 –       endif
970 –    endif
971 –
972 –    del = xa(j+1) - x
973 –    h = xa(j+1) - xa(j)
974 –
975 –    a = del / h
976 –    b = 1.0E0_DP - a
977 –    c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
978 –    d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
979 –
980 –    y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
981 –
982 –    dy = (ya(j+1)-ya(j))/h &
983 –         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
984 –         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
985 –
986 –
987 –    d2y = a*yppa(j) + b*yppa(j+1)
988 –
989 –
990 –  end subroutine eam_splint
991 –
752  
993 –  subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
753  
754  
755 <    ! yp1 and ypn are the first derivatives of y at the two endpoints
997 <    ! if boundary is 'L' the lower derivative is used
998 <    ! if boundary is 'U' the upper derivative is used
999 <    ! if boundary is 'B' then both derivatives are used
1000 <    ! if boundary is anything else, then both derivatives are assumed to be 0
1001 <
1002 <    integer :: nx, i, k, max_array_size
1003 <
1004 <    real( kind = DP ), dimension(:)        :: xa
1005 <    real( kind = DP ), dimension(:)        :: ya
1006 <    real( kind = DP ), dimension(:)        :: yppa
1007 <    real( kind = DP ), dimension(size(xa)) :: u
1008 <    real( kind = DP ) :: yp1,ypn,un,qn,sig,p
1009 <    character(len=*) :: boundary
1010 <
1011 <    ! make sure the sizes match
1012 <    if ((nx /= size(xa)) .or. (nx /= size(ya))) then
1013 <       call handleWarning("EAM_SPLINE","Array size mismatch")
1014 <    end if
1015 <
1016 <    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
1017 <         (boundary.eq.'b').or.(boundary.eq.'B')) then
1018 <       yppa(1) = -0.5E0_DP
1019 <       u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-&
1020 <            ya(1))/(xa(2)-xa(1))-yp1)
1021 <    else
1022 <       yppa(1) = 0.0E0_DP
1023 <       u(1)  = 0.0E0_DP
1024 <    endif
1025 <
1026 <    do i = 2, nx - 1
1027 <       sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1028 <       p = sig * yppa(i-1) + 2.0E0_DP
1029 <       yppa(i) = (sig - 1.0E0_DP) / p
1030 <       u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - &
1031 <            (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1032 <            (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1033 <    enddo
1034 <
1035 <    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1036 <         (boundary.eq.'b').or.(boundary.eq.'B')) then
1037 <       qn = 0.5E0_DP
1038 <       un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* &
1039 <            (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1)))
1040 <    else
1041 <       qn = 0.0E0_DP
1042 <       un = 0.0E0_DP
1043 <    endif
1044 <
1045 <    yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1046 <
1047 <    do k = nx-1, 1, -1
1048 <       yppa(k)=yppa(k)*yppa(k+1)+u(k)
1049 <    enddo
1050 <
1051 <  end subroutine eam_spline
1052 <
1053 < end module eam
755 > end module suttonchen

Diff Legend

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