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 1150 by gezelter, Fri May 7 21:35:05 2004 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 :: rho_tmp
70 >  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col
71 >  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
72   #endif
73  
74    type, private :: EAMTypeList
# Line 78 | Line 79 | module eam
79    end type EAMTypeList
80  
81  
82 <  type (eamTypeList) :: EAMList
82 >  type (eamTypeList), save :: EAMList
83  
84    !! standard eam stuff  
85  
86  
87    public :: init_EAM_FF
88 < !  public :: EAM_new_rcut
88 >  public :: setCutoffEAM
89    public :: do_eam_pair
90    public :: newEAMtype
91    public :: calc_eam_prepair_rho
92    public :: calc_eam_preforce_Frho
93 <  
93 >  public :: clean_EAM
94  
95   contains
96  
# Line 114 | Line 115 | contains
115      integer                                :: alloc_stat
116      integer                                :: current
117      integer,pointer                        :: Matchlist(:) => null()
118 <    type (EAMtype), pointer                :: makeEamtype => null()
118 >
119      status = 0
120  
121 +
122      !! Assume that atypes has already been set and get the total number of types in atypes
123      !! Also assume that every member of atypes is a EAM model.
124    
# Line 132 | Line 134 | contains
134      current = EAMList%currentAddition
135      
136  
137 <    call allocate_EAMType(eam_nrho,eam_nr,makeEamtype,stat=alloc_stat)
137 >    call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
138      if (alloc_stat /= 0) then
139         status = -1
140         return
141      end if
140 –    makeEamtype => EAMList%EAMParams(current)
142  
143      ! this is a possible bug, we assume a correspondence between the vector atypes and
144      ! EAMAtypes
# Line 147 | Line 148 | contains
148      EAMList%EAMParams(current)%eam_nrho     = eam_nrho
149      EAMList%EAMParams(current)%eam_drho     = eam_drho
150      EAMList%EAMParams(current)%eam_nr       = eam_nr
151 +    EAMList%EAMParams(current)%eam_dr       = eam_dr
152      EAMList%EAMParams(current)%eam_rcut     = rcut
153      EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
154      EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
# Line 164 | Line 166 | contains
166      integer :: number_r, number_rho
167  
168  
169 +    status = 0
170 +    if (EAMList%currentAddition == 0) then
171 +       call handleError("init_EAM_FF","No members in EAMList")
172 +       status = -1
173 +       return
174 +    end if
175  
176 +
177         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
178  
179 <          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
179 > ! Build array of r values
180  
181 +          do j = 1,EAMList%EAMParams(i)%eam_nr
182 +             EAMList%EAMParams(i)%eam_rvals(j) = &
183 +                  real(j-1,kind=dp)* &
184 +                  EAMList%EAMParams(i)%eam_dr
185 +              end do
186 + ! Build array of rho values
187 +          do j = 1,EAMList%EAMParams(i)%eam_nrho
188 +             EAMList%EAMParams(i)%eam_rhovals(j) = &
189 +                  real(j-1,kind=dp)* &
190 +                  EAMList%EAMParams(i)%eam_drho
191 +          end do
192            ! convert from eV to kcal / mol:
193            EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
194  
195            ! precompute the pair potential and get it into kcal / mol:
196 <          EAMList%EAMParams(i)%eam_phi_r = 0.0E0_DP
196 >          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
197            do j = 2, EAMList%EAMParams(i)%eam_nr
198               EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
199               EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
200            enddo
201         end do
202 +      
203  
204         do i = 1,  EAMList%currentAddition
205            number_r   = EAMList%EAMParams(i)%eam_nr
# Line 207 | Line 222 | contains
222                 EAMList%EAMParams(i)%eam_phi_r_pp, &
223                 0.0E0_DP, 0.0E0_DP, 'N')
224         enddo
225 <      
226 <       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
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
222 –
223 –
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 235 | Line 250 | contains
250    subroutine allocateEAM(status)
251      integer, intent(out) :: status
252  
238 –    integer :: nlocal
253   #ifdef IS_MPI
254      integer :: nrow
255      integer :: ncol
# Line 243 | Line 257 | contains
257      integer :: alloc_stat
258  
259  
260 <    nlocal = getNlocal()
247 <
260 >    status = 0
261   #ifdef IS_MPI
262      nrow = getNrow(plan_row)
263      ncol = getNcol(plan_col)
# Line 279 | Line 292 | contains
292      
293   #ifdef IS_MPI
294  
295 +    if (allocated(rho_tmp)) deallocate(rho_tmp)
296 +    allocate(rho_tmp(nlocal),stat=alloc_stat)
297 +    if (alloc_stat /= 0) then
298 +       status = -1
299 +       return
300 +    end if
301 +
302 +
303      if (allocated(frho_row)) deallocate(frho_row)
304      allocate(frho_row(nrow),stat=alloc_stat)
305      if (alloc_stat /= 0) then
# Line 336 | Line 357 | contains
357  
358    end subroutine allocateEAM
359  
360 + !! C sets rcut to be the largest cutoff of any atype
361 + !! present in this simulation. Doesn't include all atypes
362 + !! sim knows about, just those in the simulation.
363 +  subroutine setCutoffEAM(rcut, status)
364 +    real(kind=dp) :: rcut
365 +    integer :: status
366 +    status = 0
367  
368 <  subroutine clean_EAM()
368 >    EAM_rcut = rcut
369  
370 < ! clean non-IS_MPI first
370 >  end subroutine setCutoffEAM
371 >
372 >
373 >
374 >  subroutine clean_EAM()
375 >  
376 >   ! clean non-IS_MPI first
377      frho = 0.0_dp
378      rho  = 0.0_dp
379      dfrhodrho = 0.0_dp
# Line 349 | Line 383 | contains
383      frho_col = 0.0_dp
384      rho_row  = 0.0_dp
385      rho_col  = 0.0_dp
386 +    rho_tmp  = 0.0_dp
387      dfrhodrho_row = 0.0_dp
388      dfrhodrho_col = 0.0_dp
389   #endif
# Line 359 | Line 394 | contains
394    subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
395      integer, intent(in)          :: eam_n_rho
396      integer, intent(in)          :: eam_n_r
397 <    type (EAMType), pointer      :: thisEAMType
397 >    type (EAMType)               :: thisEAMType
398      integer, optional   :: stat
399      integer             :: alloc_stat
400  
# Line 458 | Line 493 | contains
493      integer :: myid_atom2
494  
495   ! check to see if we need to be cleaned at the start of a force loop
496 <    if (cleanme) call clean_EAM
462 <    cleanme = .false.
463 <    
496 >      
497  
498 +
499 +
500   #ifdef IS_MPI
501      myid_atom1 = atid_Row(atom1)
502      myid_atom2 = atid_Col(atom2)
# Line 471 | Line 506 | contains
506   #endif
507  
508      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
474 –
509  
510  
511 +
512         call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
513              EAMList%EAMParams(myid_atom1)%eam_rvals, &
514              EAMList%EAMParams(myid_atom1)%eam_rho_r, &
# Line 487 | Line 522 | contains
522   #else
523         rho(atom2) = rho(atom2) + rho_i_at_j
524   #endif
525 + !       write(*,*) atom1,atom2,r,rho_i_at_j
526         endif
527  
528         if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
# Line 496 | Line 532 | contains
532                 EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
533                 r, rho_j_at_i,drho,d2rho)
534  
535 <
535 >    
536        
537        
538   #ifdef  IS_MPI
# Line 506 | Line 542 | contains
542   #endif
543         endif
544  
545 +
546 +
547 +
548 +
549 +
550    end subroutine calc_eam_prepair_rho
551  
552  
# Line 521 | Line 562 | contains
562      integer :: atype1
563      integer :: me
564      integer :: n_rho_points
524 –    ! reset clean forces to be true at top of calc rho.
525 –    cleanme = .true.
565  
566 < !! Scatter the electron density from  pre-pair calculation back to local atoms
566 >  
567 >    cleanme = .true.
568 >    !! Scatter the electron density from  pre-pair calculation back to local atoms
569   #ifdef IS_MPI
570      call scatter(rho_row,rho,plan_row,eam_err)
571      if (eam_err /= 0 ) then
572        write(errMsg,*) " Error scattering rho_row into rho"
573        call handleError(RoutineName,errMesg)
574     endif      
575 <    call scatter(rho_col,rho,plan_col,eam_err)
575 >    call scatter(rho_col,rho_tmp,plan_col,eam_err)
576      if (eam_err /= 0 ) then
577        write(errMsg,*) " Error scattering rho_col into rho"
578        call handleError(RoutineName,errMesg)
579     endif
580 +
581 +      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
582   #endif
583  
584  
585 +
586   !! Calculate F(rho) and derivative
587      do atom = 1, nlocal
588         me = atid(atom)
# Line 562 | Line 606 | contains
606         end if
607  
608  
609 <       frho(i) = u
610 <       dfrhodrho(i) = u1
611 <       d2frhodrhodrho(i) = u2
609 >       frho(atom) = u
610 >       dfrhodrho(atom) = u1
611 >       d2frhodrhodrho(atom) = u2
612         pot = pot + u
613 +
614      enddo
615  
616 +  
617  
572 –
618   #ifdef IS_MPI
619      !! communicate f(rho) and derivatives back into row and column arrays
620      call gather(frho,frho_row,plan_row, eam_err)
# Line 599 | Line 644 | contains
644      endif
645   #endif
646  
647 +  
648    end subroutine calc_eam_preforce_Frho
649  
650  
651  
652  
653    !! Does EAM pairwise Force calculation.  
654 <  subroutine do_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress)
654 >  subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, pot, f, &
655 >       do_pot, do_stress)
656      !Arguments    
657      integer, intent(in) ::  atom1, atom2
658      real( kind = dp ), intent(in) :: rij, r2
659 <    real( kind = dp ) :: pot
660 <    real( kind = dp ), dimension(3,getNlocal()) :: f
659 >    real( kind = dp ) :: pot, sw, vpair
660 >    real( kind = dp ), dimension(3,nLocal) :: f
661      real( kind = dp ), intent(in), dimension(3) :: d
662      logical, intent(in) :: do_pot, do_stress
663      
# Line 631 | Line 678 | contains
678      integer  :: mytype_atom1
679      integer  :: mytype_atom2
680  
634 –
681   !Local Variables
682      
683 +   ! write(*,*) "Frho: ", Frho(atom1)
684  
638 –    
685      phab = 0.0E0_DP
686      dvpdr = 0.0E0_DP
687      d2vpdrdr = 0.0E0_DP
688 <    
643 <    
688 >
689      if (rij .lt. EAM_rcut) then
690   #ifdef IS_MPI
691   !!!!! FIX ME
# Line 727 | Line 772 | contains
772   #ifdef IS_MPI
773         dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
774              + dvpdr
775 <    
775 >
776   #else
777         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
778              + dvpdr
779 <
779 >      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
780   #endif
781  
782         fx = dudr * drdx
# Line 744 | Line 789 | contains
789            pot_Row(atom1) = pot_Row(atom1) + phab*0.5
790            pot_Col(atom2) = pot_Col(atom2) + phab*0.5
791         end if
792 +       vpair = vpair + phab
793  
794         f_Row(1,atom1) = f_Row(1,atom1) + fx
795         f_Row(2,atom1) = f_Row(2,atom1) + fy
# Line 753 | Line 799 | contains
799         f_Col(2,atom2) = f_Col(2,atom2) - fy
800         f_Col(3,atom2) = f_Col(3,atom2) - fz
801   #else
802 <       if(do_pot) pot = pot + phab
803 <      
802 >
803 >       if(do_pot) then
804 >          pot = pot + phab
805 >       end if
806 >       vpair = vpair + phab
807 >
808         f(1,atom1) = f(1,atom1) + fx
809         f(2,atom1) = f(2,atom1) + fy
810         f(3,atom1) = f(3,atom1) + fz
# Line 802 | Line 852 | contains
852   #endif
853            
854            if (molMembershipList(id1) .ne. molMembershipList(id2)) then
805 –            
855  
807 –
808 –
856               tau_Temp(1) = tau_Temp(1) - d(1) * fx
857               tau_Temp(2) = tau_Temp(2) - d(1) * fy
858               tau_Temp(3) = tau_Temp(3) - d(1) * fz

Diff Legend

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