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 650 by chuckv, Thu Jul 24 19:57:35 2003 UTC vs.
Revision 669 by chuckv, Thu Aug 7 00:47:33 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 51 | Line 51 | module eam
51  
52  
53    !! Arrays for derivatives used in force calculation
54 <  real( kind = dp), dimension(:), allocatable :: frho
55 <  real( kind = dp), dimension(:), allocatable :: rho
54 >  real( kind = dp),save, dimension(:), allocatable :: frho
55 >  real( kind = dp),save, dimension(:), allocatable :: rho
56  
57 <  real( kind = dp), dimension(:), allocatable :: dfrhodrho
58 <  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
57 >  real( kind = dp),save, dimension(:), allocatable :: dfrhodrho
58 >  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho
59  
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 84 | Line 84 | module eam
84  
85  
86    public :: init_EAM_FF
87 < !  public :: EAM_new_rcut
87 >  public :: setCutoffEAM
88    public :: do_eam_pair
89    public :: newEAMtype
90    public :: calc_eam_prepair_rho
# Line 114 | Line 114 | contains
114      integer                                :: alloc_stat
115      integer                                :: current
116      integer,pointer                        :: Matchlist(:) => null()
117 <    type (EAMtype), pointer                :: makeEamtype => null()
117 >
118      status = 0
119  
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 132 | Line 133 | contains
133      current = EAMList%currentAddition
134      
135  
136 <    call allocate_EAMType(eam_nrho,eam_nr,makeEamtype,stat=alloc_stat)
136 >    call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
137      if (alloc_stat /= 0) then
138         status = -1
139         return
140      end if
140 –    makeEamtype => EAMList%EAMParams(current)
141  
142      ! this is a possible bug, we assume a correspondence between the vector atypes and
143      ! EAMAtypes
# Line 147 | Line 147 | contains
147      EAMList%EAMParams(current)%eam_nrho     = eam_nrho
148      EAMList%EAMParams(current)%eam_drho     = eam_drho
149      EAMList%EAMParams(current)%eam_nr       = eam_nr
150 +    EAMList%EAMParams(current)%eam_dr       = eam_dr
151      EAMList%EAMParams(current)%eam_rcut     = rcut
152      EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
153      EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
# Line 163 | Line 164 | contains
164      integer :: alloc_stat
165      integer :: number_r, number_rho
166  
167 +    if (EAMList%currentAddition == 0) then
168 +       call handleError("init_EAM_FF","No members in EAMList")
169 +       status = -1
170 +       return
171 +    end if
172  
173  
174 +
175         do i = 1, EAMList%currentAddition
169 –  
170 –          EAMList%EAMParams(i)%eam_rvals(1:EAMList%EAMParams(i)%eam_nr) = &
171 –               real(EAMList%EAMParams(i)%eam_nr,kind=dp)* &
172 –               EAMList%EAMParams(i)%eam_dr
176  
177 <          EAMList%EAMParams(i)%eam_rhovals(1:EAMList%EAMParams(i)%eam_nrho) = &
175 <               real(EAMList%EAMParams(i)%eam_nrho,kind=dp)* &
176 <               EAMList%EAMParams(i)%eam_drho
177 > ! Build array of r values
178  
179 +          do j = 1,EAMList%EAMParams(i)%eam_nr
180 +             EAMList%EAMParams(i)%eam_rvals(j) = &
181 +                  real(j-1,kind=dp)* &
182 +                  EAMList%EAMParams(i)%eam_dr
183 +              end do
184 + ! Build array of rho values
185 +          do j = 1,EAMList%EAMParams(i)%eam_nrho
186 +             EAMList%EAMParams(i)%eam_rhovals(j) = &
187 +                  real(j-1,kind=dp)* &
188 +                  EAMList%EAMParams(i)%eam_drho
189 +          end do
190            ! convert from eV to kcal / mol:
191            EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
192  
193            ! precompute the pair potential and get it into kcal / mol:
194 <          EAMList%EAMParams(i)%eam_phi_r = 0.0E0_DP
194 >          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
195            do j = 2, EAMList%EAMParams(i)%eam_nr
196               EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
197               EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
198            enddo
199         end do
200 +      
201  
202         do i = 1,  EAMList%currentAddition
203            number_r   = EAMList%EAMParams(i)%eam_nr
# Line 207 | Line 220 | contains
220                 EAMList%EAMParams(i)%eam_phi_r_pp, &
221                 0.0E0_DP, 0.0E0_DP, 'N')
222         enddo
223 <      
224 <       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
223 >
224 >
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
# Line 336 | Line 350 | contains
350  
351    end subroutine allocateEAM
352  
353 + !! C sets rcut to be the largest cutoff of any atype
354 + !! present in this simulation. Doesn't include all atypes
355 + !! sim knows about, just those in the simulation.
356 +  subroutine setCutoffEAM(rcut, status)
357 +    real(kind=dp) :: rcut
358 +    integer :: status
359 +    status = 0
360  
361 <  subroutine clean_EAM()
361 >    EAM_rcut = rcut
362  
363 < ! clean non-IS_MPI first
363 >  end subroutine setCutoffEAM
364 >
365 >
366 >
367 >  subroutine clean_EAM()
368 >   ! clean non-IS_MPI first
369      frho = 0.0_dp
370      rho  = 0.0_dp
371      dfrhodrho = 0.0_dp
# Line 359 | Line 385 | contains
385    subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
386      integer, intent(in)          :: eam_n_rho
387      integer, intent(in)          :: eam_n_r
388 <    type (EAMType), pointer      :: thisEAMType
388 >    type (EAMType)               :: thisEAMType
389      integer, optional   :: stat
390      integer             :: alloc_stat
391  
# Line 458 | Line 484 | contains
484      integer :: myid_atom2
485  
486   ! check to see if we need to be cleaned at the start of a force loop
461 –    if (cleanme) call clean_EAM
462 –    cleanme = .false.
487      
488 +    if (cleanme) then
489 +       call clean_EAM
490 +       cleanme = .false.
491 +    end if
492 +      
493  
494 +    
495 +
496   #ifdef IS_MPI
497      myid_atom1 = atid_Row(atom1)
498      myid_atom2 = atid_Col(atom2)
# Line 506 | Line 537 | contains
537   #endif
538         endif
539  
540 +    
541 +
542    end subroutine calc_eam_prepair_rho
543  
544  
# Line 523 | Line 556 | contains
556      integer :: n_rho_points
557      ! reset clean forces to be true at top of calc rho.
558      cleanme = .true.
559 <
559 >  
560   !! Scatter the electron density from  pre-pair calculation back to local atoms
561   #ifdef IS_MPI
562      call scatter(rho_row,rho,plan_row,eam_err)
# Line 568 | Line 601 | contains
601         pot = pot + u
602      enddo
603  
604 +  
605  
572 –
606   #ifdef IS_MPI
607      !! communicate f(rho) and derivatives back into row and column arrays
608      call gather(frho,frho_row,plan_row, eam_err)
# Line 635 | Line 668 | contains
668   !Local Variables
669      
670  
671 <    
671 >
672      phab = 0.0E0_DP
673      dvpdr = 0.0E0_DP
674      d2vpdrdr = 0.0E0_DP
675 <    
643 <    
675 >
676      if (rij .lt. EAM_rcut) then
677   #ifdef IS_MPI
678   !!!!! FIX ME
# Line 727 | Line 759 | contains
759   #ifdef IS_MPI
760         dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
761              + dvpdr
762 <    
762 >
763   #else
764         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
765              + dvpdr
766 <
766 >      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
767   #endif
768  
769         fx = dudr * drdx
# Line 753 | Line 785 | contains
785         f_Col(2,atom2) = f_Col(2,atom2) - fy
786         f_Col(3,atom2) = f_Col(3,atom2) - fz
787   #else
788 <       if(do_pot) pot = pot + phab
789 <      
788 >
789 >       if(do_pot) then
790 >          pot = pot + phab
791 >       end if
792 >
793         f(1,atom1) = f(1,atom1) + fx
794         f(2,atom1) = f(2,atom1) + fy
795         f(3,atom1) = f(3,atom1) + fz

Diff Legend

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