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 648 by chuckv, Wed Jul 23 22:13:59 2003 UTC vs.
Revision 801 by chuckv, Sat Oct 4 18:46:12 2003 UTC

# Line 4 | Line 4 | module eam
4    use force_globals
5    use status
6    use atype_module
7 +  use Vector_class
8   #ifdef IS_MPI
9    use mpiSimulation
10   #endif
# Line 14 | 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 59 | 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 77 | 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  
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
91    public :: calc_eam_preforce_Frho
92 <  
92 >  public :: clean_EAM
93  
94   contains
95  
# Line 113 | 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 131 | 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
139 –    makeEamtype => EAMList%EAMParams(current)
141  
142      ! this is a possible bug, we assume a correspondence between the vector atypes and
143      ! EAMAtypes
# Line 146 | 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 165 | contains
165      integer :: number_r, number_rho
166  
167  
168 +    status = 0
169 +    if (EAMList%currentAddition == 0) then
170 +       call handleError("init_EAM_FF","No members in EAMList")
171 +       status = -1
172 +       return
173 +    end if
174  
175 +
176         do i = 1, EAMList%currentAddition
168 –  
169 –          EAMList%EAMParams(i)%eam_rvals(1:EAMList%EAMParams(i)%eam_nr) = &
170 –               real(EAMList%EAMParams(i)%eam_nr,kind=dp)* &
171 –               EAMList%EAMParams(i)%eam_dr
177  
178 <          EAMList%EAMParams(i)%eam_rhovals(1:EAMList%EAMParams(i)%eam_nrho) = &
174 <               real(EAMList%EAMParams(i)%eam_nrho,kind=dp)* &
175 <               EAMList%EAMParams(i)%eam_drho
178 > ! Build array of r values
179  
180 +          do j = 1,EAMList%EAMParams(i)%eam_nr
181 +             EAMList%EAMParams(i)%eam_rvals(j) = &
182 +                  real(j-1,kind=dp)* &
183 +                  EAMList%EAMParams(i)%eam_dr
184 +              end do
185 + ! Build array of rho values
186 +          do j = 1,EAMList%EAMParams(i)%eam_nrho
187 +             EAMList%EAMParams(i)%eam_rhovals(j) = &
188 +                  real(j-1,kind=dp)* &
189 +                  EAMList%EAMParams(i)%eam_drho
190 +          end do
191            ! convert from eV to kcal / mol:
192            EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
193  
194            ! precompute the pair potential and get it into kcal / mol:
195 <          EAMList%EAMParams(i)%eam_phi_r = 0.0E0_DP
195 >          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
196            do j = 2, EAMList%EAMParams(i)%eam_nr
197               EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
198               EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
199            enddo
200         end do
201 +      
202  
203         do i = 1,  EAMList%currentAddition
204            number_r   = EAMList%EAMParams(i)%eam_nr
# Line 206 | Line 221 | contains
221                 EAMList%EAMParams(i)%eam_phi_r_pp, &
222                 0.0E0_DP, 0.0E0_DP, 'N')
223         enddo
224 <      
225 <       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
224 >
225 >
226 > !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
227         !! find the smallest rcut for any eam atype
228 <       do i = 2, EAMList%currentAddition
229 <          current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
230 <       end do
228 > !       do i = 2, EAMList%currentAddition
229 > !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
230 > !       end do
231  
232 <       EAM_rcut = current_rcut_max
233 <       EAM_rcut_orig = current_rcut_max
232 > !       EAM_rcut = current_rcut_max
233 > !       EAM_rcut_orig = current_rcut_max
234   !       do i = 1, EAMList%currentAddition
235   !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
236   !       end do
221 –
222 –
237         !! Allocate arrays for force calculation
238 <          call allocateEAM(alloc_stat)
239 <          if (alloc_stat /= 0 ) then
240 <             status = -1
241 <             return
242 <          endif
243 <
238 >      
239 >       call allocateEAM(alloc_stat)
240 >       if (alloc_stat /= 0 ) then
241 >          write(*,*) "allocateEAM failed"
242 >          status = -1
243 >          return
244 >       endif
245 >    
246    end subroutine init_EAM_FF
247  
248   !! routine checks to see if array is allocated, deallocates array if allocated
# Line 243 | Line 259 | contains
259  
260  
261      nlocal = getNlocal()
262 <
262 >    status = 0
263   #ifdef IS_MPI
264      nrow = getNrow(plan_row)
265      ncol = getNcol(plan_col)
# Line 335 | Line 351 | contains
351  
352    end subroutine allocateEAM
353  
354 + !! C sets rcut to be the largest cutoff of any atype
355 + !! present in this simulation. Doesn't include all atypes
356 + !! sim knows about, just those in the simulation.
357 +  subroutine setCutoffEAM(rcut, status)
358 +    real(kind=dp) :: rcut
359 +    integer :: status
360 +    status = 0
361  
362 <  subroutine clean_EAM()
362 >    EAM_rcut = rcut
363  
364 < ! clean non-IS_MPI first
364 >  end subroutine setCutoffEAM
365 >
366 >
367 >
368 >  subroutine clean_EAM()
369 >  
370 >   ! clean non-IS_MPI first
371      frho = 0.0_dp
372      rho  = 0.0_dp
373      dfrhodrho = 0.0_dp
# Line 358 | Line 387 | contains
387    subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
388      integer, intent(in)          :: eam_n_rho
389      integer, intent(in)          :: eam_n_r
390 <    type (EAMType), pointer      :: thisEAMType
390 >    type (EAMType)               :: thisEAMType
391      integer, optional   :: stat
392      integer             :: alloc_stat
393  
# Line 457 | Line 486 | contains
486      integer :: myid_atom2
487  
488   ! check to see if we need to be cleaned at the start of a force loop
489 <    if (cleanme) call clean_EAM
461 <    cleanme = .false.
462 <    
489 >      
490  
491 +
492 +
493   #ifdef IS_MPI
494      myid_atom1 = atid_Row(atom1)
495      myid_atom2 = atid_Col(atom2)
# Line 470 | Line 499 | contains
499   #endif
500  
501      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
473 –
502  
503  
504 +
505         call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
506              EAMList%EAMParams(myid_atom1)%eam_rvals, &
507              EAMList%EAMParams(myid_atom1)%eam_rho_r, &
# Line 486 | Line 515 | contains
515   #else
516         rho(atom2) = rho(atom2) + rho_i_at_j
517   #endif
518 + !       write(*,*) atom1,atom2,r,rho_i_at_j
519         endif
520  
521         if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
# Line 505 | Line 535 | contains
535   #endif
536         endif
537  
538 +
539    end subroutine calc_eam_prepair_rho
540  
541  
# Line 520 | Line 551 | contains
551      integer :: atype1
552      integer :: me
553      integer :: n_rho_points
523 –    ! reset clean forces to be true at top of calc rho.
524 –    cleanme = .true.
554  
555 < !! Scatter the electron density from  pre-pair calculation back to local atoms
555 >  
556 >    cleanme = .true.
557 >    !! Scatter the electron density from  pre-pair calculation back to local atoms
558   #ifdef IS_MPI
559      call scatter(rho_row,rho,plan_row,eam_err)
560      if (eam_err /= 0 ) then
# Line 561 | Line 592 | contains
592         end if
593  
594  
595 <       frho(i) = u
596 <       dfrhodrho(i) = u1
597 <       d2frhodrhodrho(i) = u2
595 >       frho(atom) = u
596 >       dfrhodrho(atom) = u1
597 >       d2frhodrhodrho(atom) = u2
598         pot = pot + u
599 +
600      enddo
601  
602 +  
603  
571 –
604   #ifdef IS_MPI
605      !! communicate f(rho) and derivatives back into row and column arrays
606      call gather(frho,frho_row,plan_row, eam_err)
# Line 598 | Line 630 | contains
630      endif
631   #endif
632  
633 +  
634    end subroutine calc_eam_preforce_Frho
635  
636  
# Line 630 | Line 663 | contains
663      integer  :: mytype_atom1
664      integer  :: mytype_atom2
665  
633 –
666   !Local Variables
667      
668 +   ! write(*,*) "Frho: ", Frho(atom1)
669  
637 –    
670      phab = 0.0E0_DP
671      dvpdr = 0.0E0_DP
672      d2vpdrdr = 0.0E0_DP
673 <    
642 <    
673 >
674      if (rij .lt. EAM_rcut) then
675   #ifdef IS_MPI
676   !!!!! FIX ME
# Line 726 | Line 757 | contains
757   #ifdef IS_MPI
758         dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
759              + dvpdr
760 <    
760 >
761   #else
762         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
763              + dvpdr
764 <
764 >      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
765   #endif
766  
767         fx = dudr * drdx
# Line 752 | Line 783 | contains
783         f_Col(2,atom2) = f_Col(2,atom2) - fy
784         f_Col(3,atom2) = f_Col(3,atom2) - fz
785   #else
786 <       if(do_pot) pot = pot + phab
787 <      
786 >
787 >       if(do_pot) then
788 >          pot = pot + phab
789 >       end if
790 >
791         f(1,atom1) = f(1,atom1) + fx
792         f(2,atom1) = f(2,atom1) + fy
793         f(3,atom1) = f(3,atom1) + fz
# Line 801 | Line 835 | contains
835   #endif
836            
837            if (molMembershipList(id1) .ne. molMembershipList(id2)) then
804 –            
838  
806 –
807 –
839               tau_Temp(1) = tau_Temp(1) - d(1) * fx
840               tau_Temp(2) = tau_Temp(2) - d(1) * fy
841               tau_Temp(3) = tau_Temp(3) - d(1) * fz

Diff Legend

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