ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_eam.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/calc_eam.F90 (file contents):
Revision 657 by chuckv, Wed Jul 30 21:17:01 2003 UTC vs.
Revision 747 by gezelter, Fri Sep 5 21:28:52 2003 UTC

# Line 15 | Line 15 | module eam
15    logical, save :: EAM_FF_initialized = .false.
16    integer, save :: EAM_Mixing_Policy
17    real(kind = dp), save :: EAM_rcut
18 <  real(kind = dp), save :: EAM_rcut_orig
18 >  logical, save :: haveRcut = .false.
19  
20    character(len = statusMsgSize) :: errMesg
21    integer :: eam_err
# Line 60 | Line 60 | module eam
60  
61   !! Arrays for MPI storage
62   #ifdef IS_MPI
63 <  real( kind = dp), dimension(:), allocatable :: dfrhodrho_col
64 <  real( kind = dp), dimension(:), allocatable :: dfrhodrho_row
65 <  real( kind = dp), dimension(:), allocatable :: frho_row
66 <  real( kind = dp), dimension(:), allocatable :: frho_col
67 <  real( kind = dp), dimension(:), allocatable :: rho_row
68 <  real( kind = dp), dimension(:), allocatable :: rho_col
69 <  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_col
70 <  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_row
63 >  real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
64 >  real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row
65 >  real( kind = dp),save, dimension(:), allocatable :: frho_row
66 >  real( kind = dp),save, dimension(:), allocatable :: frho_col
67 >  real( kind = dp),save, dimension(:), allocatable :: rho_row
68 >  real( kind = dp),save, dimension(:), allocatable :: rho_col
69 >  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col
70 >  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
71   #endif
72  
73    type, private :: EAMTypeList
# Line 78 | Line 78 | module eam
78    end type EAMTypeList
79  
80  
81 <  type (eamTypeList) :: EAMList
81 >  type (eamTypeList), save :: EAMList
82  
83    !! standard eam stuff  
84  
# Line 89 | Line 89 | module eam
89    public :: newEAMtype
90    public :: calc_eam_prepair_rho
91    public :: calc_eam_preforce_Frho
92 <  
92 >  public :: clean_EAM
93  
94   contains
95  
# Line 117 | Line 117 | contains
117  
118      status = 0
119  
120 <    write(*,*) "Adding new eamtype: ",eam_ident
120 >
121      !! Assume that atypes has already been set and get the total number of types in atypes
122      !! Also assume that every member of atypes is a EAM model.
123    
# Line 222 | Line 222 | contains
222         enddo
223  
224  
225 <       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
225 > !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
226         !! find the smallest rcut for any eam atype
227 <       do i = 2, EAMList%currentAddition
228 <          current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
229 <       end do
227 > !       do i = 2, EAMList%currentAddition
228 > !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
229 > !       end do
230  
231 <       EAM_rcut = current_rcut_max
232 <       EAM_rcut_orig = current_rcut_max
231 > !       EAM_rcut = current_rcut_max
232 > !       EAM_rcut_orig = current_rcut_max
233   !       do i = 1, EAMList%currentAddition
234   !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
235   !       end do
236 –
237 –
236         !! Allocate arrays for force calculation
237            call allocateEAM(alloc_stat)
238            if (alloc_stat /= 0 ) then
# Line 350 | Line 348 | contains
348  
349    end subroutine allocateEAM
350  
351 + !! C sets rcut to be the largest cutoff of any atype
352 + !! present in this simulation. Doesn't include all atypes
353 + !! sim knows about, just those in the simulation.
354    subroutine setCutoffEAM(rcut, status)
355      real(kind=dp) :: rcut
356      integer :: status
357 +    status = 0
358  
359 <    if (rcut < EAM_rcut) then
358 <       EAM_rcut = rcut
359 <    endif
359 >    EAM_rcut = rcut
360  
361 –
361    end subroutine setCutoffEAM
362  
363  
364  
365    subroutine clean_EAM()
366 <
367 < ! clean non-IS_MPI first
366 >  
367 >   ! clean non-IS_MPI first
368      frho = 0.0_dp
369      rho  = 0.0_dp
370      dfrhodrho = 0.0_dp
# Line 484 | Line 483 | contains
483      integer :: myid_atom2
484  
485   ! check to see if we need to be cleaned at the start of a force loop
486 <    if (cleanme) call clean_EAM
488 <    cleanme = .false.
489 <    
486 >      
487  
488 +
489 +
490   #ifdef IS_MPI
491      myid_atom1 = atid_Row(atom1)
492      myid_atom2 = atid_Col(atom2)
# Line 497 | Line 496 | contains
496   #endif
497  
498      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
500 –
499  
500  
501 +
502         call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
503              EAMList%EAMParams(myid_atom1)%eam_rvals, &
504              EAMList%EAMParams(myid_atom1)%eam_rho_r, &
# Line 513 | Line 512 | contains
512   #else
513         rho(atom2) = rho(atom2) + rho_i_at_j
514   #endif
515 + !       write(*,*) atom1,atom2,r,rho_i_at_j
516         endif
517  
518         if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
# Line 532 | Line 532 | contains
532   #endif
533         endif
534  
535 +
536    end subroutine calc_eam_prepair_rho
537  
538  
# Line 547 | Line 548 | contains
548      integer :: atype1
549      integer :: me
550      integer :: n_rho_points
550 –    ! reset clean forces to be true at top of calc rho.
551 –    cleanme = .true.
551  
552 < !! Scatter the electron density from  pre-pair calculation back to local atoms
552 >  
553 >    cleanme = .true.
554 >    !! Scatter the electron density from  pre-pair calculation back to local atoms
555   #ifdef IS_MPI
556      call scatter(rho_row,rho,plan_row,eam_err)
557      if (eam_err /= 0 ) then
# Line 588 | Line 589 | contains
589         end if
590  
591  
592 <       frho(i) = u
593 <       dfrhodrho(i) = u1
594 <       d2frhodrhodrho(i) = u2
592 >       frho(atom) = u
593 >       dfrhodrho(atom) = u1
594 >       d2frhodrhodrho(atom) = u2
595         pot = pot + u
595 –    enddo
596  
597 +    enddo
598  
599 +  
600  
601   #ifdef IS_MPI
602      !! communicate f(rho) and derivatives back into row and column arrays
# Line 625 | Line 627 | contains
627      endif
628   #endif
629  
630 +  
631    end subroutine calc_eam_preforce_Frho
632  
633  
# Line 657 | Line 660 | contains
660      integer  :: mytype_atom1
661      integer  :: mytype_atom2
662  
660 –
663   !Local Variables
664      
665 +   ! write(*,*) "Frho: ", Frho(atom1)
666  
664 –    
667      phab = 0.0E0_DP
668      dvpdr = 0.0E0_DP
669      d2vpdrdr = 0.0E0_DP
670 <    
669 <    
670 >
671      if (rij .lt. EAM_rcut) then
672   #ifdef IS_MPI
673   !!!!! FIX ME
# Line 753 | Line 754 | contains
754   #ifdef IS_MPI
755         dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
756              + dvpdr
757 <    
757 >
758   #else
759         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
760              + dvpdr
761 <
761 >      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
762   #endif
763  
764         fx = dudr * drdx
# Line 779 | Line 780 | contains
780         f_Col(2,atom2) = f_Col(2,atom2) - fy
781         f_Col(3,atom2) = f_Col(3,atom2) - fz
782   #else
783 <       if(do_pot) pot = pot + phab
784 <      
783 >
784 >       if(do_pot) then
785 >          pot = pot + phab
786 >       end if
787 >
788         f(1,atom1) = f(1,atom1) + fx
789         f(2,atom1) = f(2,atom1) + fy
790         f(3,atom1) = f(3,atom1) + fz
# Line 828 | Line 832 | contains
832   #endif
833            
834            if (molMembershipList(id1) .ne. molMembershipList(id2)) then
831 –            
835  
833 –
834 –
836               tau_Temp(1) = tau_Temp(1) - d(1) * fx
837               tau_Temp(2) = tau_Temp(2) - d(1) * fy
838               tau_Temp(3) = tau_Temp(3) - d(1) * fz

Diff Legend

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