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

Comparing trunk/OOPSE-2.0/src/UseTheForce/DarkSide/suttonchen.F90 (file contents):
Revision 2412 by chuckv, Wed Nov 2 23:50:56 2005 UTC vs.
Revision 2413 by chuckv, Thu Nov 3 23:06:00 2005 UTC

# Line 66 | Line 66 | module suttonchen
66    logical, save,:: haveMixingMap = .false.
67  
68    character(len = statusMsgSize) :: errMesg
69 <  integer :: eam_err
69 >  integer :: sc_err
70  
71    character(len = 200) :: errMsg
72    character(len=*), parameter :: RoutineName =  "Sutton-Chen MODULE"
# Line 174 | Line 174 | contains
174  
175  
176      ! check to see if this is the first time into
177 <    if (.not.associated(EAMList%EAMParams)) then
177 >    if (.not.associated(SCTypeList%EAMParams)) then
178         call getMatchingElementList(atypes, "is_SuttonChen", .true., nSCtypes, MatchList)
179         SCTypeList%nSCtypes = nSCtypes
180         allocate(SCTypeList%SCTypes(nSCTypes))
# Line 215 | Line 215 | contains
215      integer :: eamID
216      real(kind=dp) :: cutValue
217      
218 <    eamID = EAMList%atidToEAMType(atomID)
219 <    cutValue = EAMList%EAMParams(eamID)%eam_rcut
218 >    eamID = SCTypeList%atidToEAMType(atomID)
219 >    cutValue = SCTypeList%EAMParams(eamID)%eam_rcut
220    end function getSCCut
221  
222  
# Line 432 | Line 432 | contains
432  
433      ! we don't use the derivatives, dummy variables
434      real( kind = dp) :: drho,d2rho
435 <    integer :: eam_err
435 >    integer :: sc_err
436      
437      integer :: atid1,atid2 ! Global atid    
438      integer :: myid_atom1 ! EAM atid
# Line 469 | Line 469 | contains
469         rho(atom1) = rho(atom1) + rho_j_at_i
470   #endif
471  
472 –
473 –
472    end subroutine calc_sc_prepair_rho
473  
474  
475 < !! Calculate the functional F(rho) for all local atoms
475 > !! Calculate the rho_a for all local atoms
476    subroutine calc_eam_preforce_Frho(nlocal,pot)
477      integer :: nlocal
478      real(kind=dp) :: pot
# Line 489 | Line 487 | contains
487      cleanme = .true.
488      !! Scatter the electron density from  pre-pair calculation back to local atoms
489   #ifdef IS_MPI
490 <    call scatter(rho_row,rho,plan_atom_row,eam_err)
491 <    if (eam_err /= 0 ) then
490 >    call scatter(rho_row,rho,plan_atom_row,sc_err)
491 >    if (sc_err /= 0 ) then
492         write(errMsg,*) " Error scattering rho_row into rho"
493         call handleError(RoutineName,errMesg)
494      endif
495 <    call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
496 <    if (eam_err /= 0 ) then
495 >    call scatter(rho_col,rho_tmp,plan_atom_col,sc_err)
496 >    if (sc_err /= 0 ) then
497         write(errMsg,*) " Error scattering rho_col into rho"
498         call handleError(RoutineName,errMesg)
499  
# Line 518 | Line 516 | contains
516  
517      #ifdef IS_MPI
518      !! communicate f(rho) and derivatives back into row and column arrays
519 <    call gather(frho,frho_row,plan_atom_row, eam_err)
520 <    if (eam_err /=  0) then
519 >    call gather(frho,frho_row,plan_atom_row, sc_err)
520 >    if (sc_err /=  0) then
521         call handleError("cal_eam_forces()","MPI gather frho_row failure")
522      endif
523 <    call gather(dfrhodrho,dfrhodrho_row,plan_atom_row, eam_err)
524 <    if (eam_err /=  0) then
523 >    call gather(dfrhodrho,dfrhodrho_row,plan_atom_row, sc_err)
524 >    if (sc_err /=  0) then
525         call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure")
526      endif
527 <    call gather(frho,frho_col,plan_atom_col, eam_err)
528 <    if (eam_err /=  0) then
527 >    call gather(frho,frho_col,plan_atom_col, sc_err)
528 >    if (sc_err /=  0) then
529         call handleError("cal_eam_forces()","MPI gather frho_col failure")
530      endif
531 <    call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, eam_err)
532 <    if (eam_err /=  0) then
531 >    call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, sc_err)
532 >    if (sc_err /=  0) then
533         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
534      endif
535  
# Line 564 | Line 562 | contains
562      real( kind = dp ) :: rha,drha,d2rha, dpha
563      real( kind = dp ) :: rhb,drhb,d2rhb, dphb
564      real( kind = dp ) :: dudr
565 <    real( kind = dp ) :: rci,rcj
565 >    real( kind = dp ) :: rcij
566      real( kind = dp ) :: drhoidr,drhojdr
567      real( kind = dp ) :: d2rhoidrdr
568      real( kind = dp ) :: d2rhojdrdr
# Line 583 | Line 581 | contains
581      dvpdr = 0.0E0_DP
582      d2vpdrdr = 0.0E0_DP
583  
586 –    if (rij .lt. EAM_rcut) then
584  
585   #ifdef IS_MPI
586         atid1 = atid_row(atom1)
# Line 593 | Line 590 | contains
590         atid2 = atid(atom2)
591   #endif
592  
593 <       mytype_atom1 = EAMList%atidToEAMType(atid1)
594 <       mytype_atom2 = EAMList%atidTOEAMType(atid2)
593 >       mytype_atom1 = SCTypeList%atidToSCType(atid1)
594 >       mytype_atom2 = SCTypeList%atidTOSCType(atid2)
595  
599 –
600 –       ! get cutoff for atom 1
601 –       rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
602 –       ! get type specific cutoff for atom 2
603 –       rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
604 –
596         drdx = d(1)/rij
597         drdy = d(2)/rij
598         drdz = d(3)/rij
599  
600 <       if (rij.lt.rci) then
601 <          call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
602 <               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
603 <               EAMList%EAMParams(mytype_atom1)%eam_rho_r, &
604 <               EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, &
605 <               rij, rha,drha,d2rha)
615 <          !! Calculate Phi(r) for atom1.
616 <          call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
617 <               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
618 <               EAMList%EAMParams(mytype_atom1)%eam_phi_r, &
619 <               EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, &
620 <               rij, pha,dpha,d2pha)
621 <       endif
622 <
623 <       if (rij.lt.rcj) then      
624 <          ! Calculate rho,drho and d2rho for atom1
625 <          call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
626 <               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
627 <               EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
628 <               EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
629 <               rij, rhb,drhb,d2rhb)
600 >                
601 >       epsilonij = MixingMap(mytype_atom1,mytype_atom2)%epsilon
602 >       aij       = MixingMap(mytype_atom1,mytype_atom2)%alpha
603 >       nij       = MixingMap(mytype_atom1,mytype_atom2)%n
604 >       mij       = MixingMap(mytype_atom1,mytype_atom2)%m
605 >       vcij      = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot              
606  
607 <          !! Calculate Phi(r) for atom2.
632 <          call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
633 <               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
634 <               EAMList%EAMParams(mytype_atom2)%eam_phi_r, &
635 <               EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, &
636 <               rij, phb,dphb,d2phb)
637 <       endif
607 >       vptmp = dij*((aij/r)**nij)
608  
639 –       if (rij.lt.rci) then
640 –          phab = phab + 0.5E0_DP*(rhb/rha)*pha
641 –          dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
642 –               pha*((drhb/rha) - (rhb*drha/rha/rha)))
643 –          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rhb/rha)*d2pha + &
644 –               2.0E0_DP*dpha*((drhb/rha) - (rhb*drha/rha/rha)) + &
645 –               pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
646 –               (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
647 –       endif
609  
610 <       if (rij.lt.rcj) then
611 <          phab = phab + 0.5E0_DP*(rha/rhb)*phb
612 <          dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
613 <               phb*((drha/rhb) - (rha*drhb/rhb/rhb)))
614 <          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rha/rhb)*d2phb + &
615 <               2.0E0_DP*dphb*((drha/rhb) - (rha*drhb/rhb/rhb)) + &
616 <               phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
617 <               (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
618 <       endif
610 >       dvpdr = -nij*vptmp/r
611 >       d2vpdrdr = -dvpdr*(nij+1.0d0)/r
612 >      
613 >       drhodr = -mij*((aij/r)**mij)/r
614 >       d2rhodrdr = -drhodr*(mij+1.0d0)/r
615 >      
616 >       dudr = drhodr*(dfrhodrho(i)+dfrhodrho(j)) &
617 >            + dvpdr
618 >
619 >
620  
659 –       drhoidr = drha
660 –       drhojdr = drhb
661 –
662 –       d2rhoidrdr = d2rha
663 –       d2rhojdrdr = d2rhb
664 –
665 –
621   #ifdef IS_MPI
622 <       dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
622 >       dudr = drhodr*(dfrhodrho_row(atom1)+dfrhodrho_col(atom2)) &
623              + dvpdr
624  
625   #else
626 <       dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
626 >       dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) &
627              + dvpdr
673 –       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
628   #endif
629  
630 +
631         fx = dudr * drdx
632         fy = dudr * drdy
633         fz = dudr * drdz

Diff Legend

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