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 669 by chuckv, Thu Aug 7 00:47:33 2003 UTC vs.
Revision 882 by chuckv, Wed Dec 17 20:13:33 2003 UTC

# Line 51 | Line 51 | module eam
51  
52  
53    !! Arrays for derivatives used in force calculation
54 <  real( kind = dp),save, dimension(:), allocatable :: frho
55 <  real( kind = dp),save, dimension(:), allocatable :: rho
54 >  real( kind = dp), dimension(:), allocatable :: frho
55 >  real( kind = dp), dimension(:), allocatable :: rho
56  
57 <  real( kind = dp),save, dimension(:), allocatable :: dfrhodrho
58 <  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho
57 >  real( kind = dp), dimension(:), allocatable :: dfrhodrho
58 >  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
59  
60  
61   !! Arrays for MPI storage
# Line 66 | Line 66 | module eam
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
# 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  
# Line 89 | Line 90 | module eam
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 164 | Line 165 | contains
165      integer :: alloc_stat
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 <
174 <
176 >
177         do i = 1, EAMList%currentAddition
178  
179   ! Build array of r values
# Line 221 | Line 223 | contains
223                 0.0E0_DP, 0.0E0_DP, 'N')
224         enddo
225  
224 –
226   !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
227         !! find the smallest rcut for any eam atype
228   !       do i = 2, EAMList%currentAddition
# Line 233 | Line 234 | contains
234   !       do i = 1, EAMList%currentAddition
235   !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
236   !       end do
236 –
237 –
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 258 | 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 293 | Line 294 | contains
294      
295   #ifdef IS_MPI
296  
297 +    if (allocated(rho_tmp)) deallocate(rho_tmp)
298 +    allocate(rho_tmp(nlocal),stat=alloc_stat)
299 +    if (alloc_stat /= 0) then
300 +       status = -1
301 +       return
302 +    end if
303 +
304 +
305      if (allocated(frho_row)) deallocate(frho_row)
306      allocate(frho_row(nrow),stat=alloc_stat)
307      if (alloc_stat /= 0) then
# Line 365 | Line 374 | contains
374  
375  
376    subroutine clean_EAM()
377 +  
378     ! clean non-IS_MPI first
379      frho = 0.0_dp
380      rho  = 0.0_dp
# Line 375 | Line 385 | contains
385      frho_col = 0.0_dp
386      rho_row  = 0.0_dp
387      rho_col  = 0.0_dp
388 +    rho_tmp  = 0.0_dp
389      dfrhodrho_row = 0.0_dp
390      dfrhodrho_col = 0.0_dp
391   #endif
# Line 484 | Line 495 | contains
495      integer :: myid_atom2
496  
497   ! check to see if we need to be cleaned at the start of a force loop
487 –    
488 –    if (cleanme) then
489 –       call clean_EAM
490 –       cleanme = .false.
491 –    end if
498        
499  
494 –    
500  
501 +
502   #ifdef IS_MPI
503      myid_atom1 = atid_Row(atom1)
504      myid_atom2 = atid_Col(atom2)
# Line 502 | Line 508 | contains
508   #endif
509  
510      if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
505 –
511  
512  
513 +
514         call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
515              EAMList%EAMParams(myid_atom1)%eam_rvals, &
516              EAMList%EAMParams(myid_atom1)%eam_rho_r, &
# Line 518 | Line 524 | contains
524   #else
525         rho(atom2) = rho(atom2) + rho_i_at_j
526   #endif
527 + !       write(*,*) atom1,atom2,r,rho_i_at_j
528         endif
529  
530         if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
# Line 527 | Line 534 | contains
534                 EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
535                 r, rho_j_at_i,drho,d2rho)
536  
537 <
537 >    
538        
539        
540   #ifdef  IS_MPI
# Line 537 | Line 544 | contains
544   #endif
545         endif
546  
540 –    
547  
548 +
549 +
550 +
551 +
552    end subroutine calc_eam_prepair_rho
553  
554  
# Line 554 | Line 564 | contains
564      integer :: atype1
565      integer :: me
566      integer :: n_rho_points
567 <    ! reset clean forces to be true at top of calc rho.
558 <    cleanme = .true.
567 >
568    
569 < !! Scatter the electron density from  pre-pair calculation back to local atoms
569 >    cleanme = .true.
570 >    !! Scatter the electron density from  pre-pair calculation back to local atoms
571   #ifdef IS_MPI
572      call scatter(rho_row,rho,plan_row,eam_err)
573      if (eam_err /= 0 ) then
574        write(errMsg,*) " Error scattering rho_row into rho"
575        call handleError(RoutineName,errMesg)
576     endif      
577 <    call scatter(rho_col,rho,plan_col,eam_err)
577 >    call scatter(rho_col,rho_tmp,plan_col,eam_err)
578      if (eam_err /= 0 ) then
579        write(errMsg,*) " Error scattering rho_col into rho"
580        call handleError(RoutineName,errMesg)
581     endif
582 +
583 +      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
584   #endif
585  
586  
587 +
588   !! Calculate F(rho) and derivative
589      do atom = 1, nlocal
590         me = atid(atom)
# Line 595 | Line 608 | contains
608         end if
609  
610  
611 <       frho(i) = u
612 <       dfrhodrho(i) = u1
613 <       d2frhodrhodrho(i) = u2
611 >       frho(atom) = u
612 >       dfrhodrho(atom) = u1
613 >       d2frhodrhodrho(atom) = u2
614         pot = pot + u
615 +
616      enddo
617  
618 <  
618 >  
619  
620   #ifdef IS_MPI
621      !! communicate f(rho) and derivatives back into row and column arrays
# Line 632 | Line 646 | contains
646      endif
647   #endif
648  
649 +  
650    end subroutine calc_eam_preforce_Frho
651  
652  
# Line 664 | Line 679 | contains
679      integer  :: mytype_atom1
680      integer  :: mytype_atom2
681  
667 –
682   !Local Variables
683      
684 <
684 >   ! write(*,*) "Frho: ", Frho(atom1)
685  
686      phab = 0.0E0_DP
687      dvpdr = 0.0E0_DP
# Line 837 | Line 851 | contains
851   #endif
852            
853            if (molMembershipList(id1) .ne. molMembershipList(id2)) then
840 –            
841 –
854  
843 –
855               tau_Temp(1) = tau_Temp(1) - d(1) * fx
856               tau_Temp(2) = tau_Temp(2) - d(1) * fy
857               tau_Temp(3) = tau_Temp(3) - d(1) * fz

Diff Legend

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