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 747 by gezelter, Fri Sep 5 21:28:52 2003 UTC vs.
Revision 1198 by tim, Thu May 27 00:48:12 2004 UTC

# 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 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 234 | Line 235 | contains
235   !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
236   !       end do
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 247 | Line 250 | contains
250    subroutine allocateEAM(status)
251      integer, intent(out) :: status
252  
250 –    integer :: nlocal
253   #ifdef IS_MPI
254 <    integer :: nrow
255 <    integer :: ncol
254 >    integer :: nAtomsInRow
255 >    integer :: nAtomsInCol
256   #endif
257      integer :: alloc_stat
258  
259  
260 <    nlocal = getNlocal()
259 <
260 >    status = 0
261   #ifdef IS_MPI
262 <    nrow = getNrow(plan_row)
263 <    ncol = getNcol(plan_col)
262 >    nAtomsInRow = getNatomsInRow(plan_atom_row)
263 >    nAtomsInCol = getNatomsInCol(plan_atom_col)
264   #endif
265  
266      if (allocated(frho)) deallocate(frho)
# Line 291 | 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)
304 >    allocate(frho_row(nAtomsInRow),stat=alloc_stat)
305      if (alloc_stat /= 0) then
306         status = -1
307         return
308      end if
309      if (allocated(rho_row)) deallocate(rho_row)
310 <    allocate(rho_row(nrow),stat=alloc_stat)
310 >    allocate(rho_row(nAtomsInRow),stat=alloc_stat)
311      if (alloc_stat /= 0) then
312         status = -1
313         return
314      end if
315      if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
316 <    allocate(dfrhodrho_row(nrow),stat=alloc_stat)
316 >    allocate(dfrhodrho_row(nAtomsInRow),stat=alloc_stat)
317      if (alloc_stat /= 0) then
318         status = -1
319         return
320      end if
321      if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
322 <    allocate(d2frhodrhodrho_row(nrow),stat=alloc_stat)
322 >    allocate(d2frhodrhodrho_row(nAtomsInRow),stat=alloc_stat)
323      if (alloc_stat /= 0) then
324         status = -1
325         return
# Line 320 | Line 329 | contains
329   ! Now do column arrays
330  
331      if (allocated(frho_col)) deallocate(frho_col)
332 <    allocate(frho_col(ncol),stat=alloc_stat)
332 >    allocate(frho_col(nAtomsInCol),stat=alloc_stat)
333      if (alloc_stat /= 0) then
334         status = -1
335         return
336      end if
337      if (allocated(rho_col)) deallocate(rho_col)
338 <    allocate(rho_col(ncol),stat=alloc_stat)
338 >    allocate(rho_col(nAtomsInCol),stat=alloc_stat)
339      if (alloc_stat /= 0) then
340         status = -1
341         return
342      end if
343      if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
344 <    allocate(dfrhodrho_col(ncol),stat=alloc_stat)
344 >    allocate(dfrhodrho_col(nAtomsInCol),stat=alloc_stat)
345      if (alloc_stat /= 0) then
346         status = -1
347         return
348      end if
349      if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
350 <    allocate(d2frhodrhodrho_col(ncol),stat=alloc_stat)
350 >    allocate(d2frhodrhodrho_col(nAtomsInCol),stat=alloc_stat)
351      if (alloc_stat /= 0) then
352         status = -1
353         return
# Line 374 | 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 522 | 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 533 | Line 543 | contains
543         endif
544  
545  
546 +
547 +
548 +
549 +
550    end subroutine calc_eam_prepair_rho
551  
552  
# Line 553 | Line 567 | contains
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)
570 >    call scatter(rho_row,rho,plan_atom_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_atom_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 600 | Line 617 | contains
617  
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)
620 >    call gather(frho,frho_row,plan_atom_row, eam_err)
621      if (eam_err /=  0) then
622         call handleError("cal_eam_forces()","MPI gather frho_row failure")
623      endif
624 <    call gather(dfrhodrho,dfrhodrho_row,plan_row, eam_err)
624 >    call gather(dfrhodrho,dfrhodrho_row,plan_atom_row, eam_err)
625      if (eam_err /=  0) then
626         call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure")
627      endif
628 <    call gather(frho,frho_col,plan_col, eam_err)
628 >    call gather(frho,frho_col,plan_atom_col, eam_err)
629      if (eam_err /=  0) then
630         call handleError("cal_eam_forces()","MPI gather frho_col failure")
631      endif
632 <    call gather(dfrhodrho,dfrhodrho_col,plan_col, eam_err)
632 >    call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, eam_err)
633      if (eam_err /=  0) then
634         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
635      endif
# Line 622 | Line 639 | contains
639  
640  
641      if (nmflag) then
642 <       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_row)
643 <       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_col)
642 >       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_atom_row)
643 >       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_atom_col)
644      endif
645   #endif
646  
# Line 634 | Line 651 | contains
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, fpair, &
655 >       pot, f, do_pot)
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
662 >    real( kind = dp ), intent(inout), dimension(3) :: fpair
663 >
664 >    logical, intent(in) :: do_pot
665      
666      real( kind = dp ) :: drdx,drdy,drdz
667      real( kind = dp ) :: d2
# Line 794 | Line 814 | contains
814         f(3,atom2) = f(3,atom2) - fz
815   #endif
816        
817 +       vpair = vpair + phab
818 + #ifdef IS_MPI
819 +       id1 = AtomRowToGlobal(atom1)
820 +       id2 = AtomColToGlobal(atom2)
821 + #else
822 +       id1 = atom1
823 +       id2 = atom2
824 + #endif
825 +      
826 +       if (molMembershipList(id1) .ne. molMembershipList(id2)) then
827 +          
828 +          fpair(1) = fpair(1) + fx
829 +          fpair(2) = fpair(2) + fy
830 +          fpair(3) = fpair(3) + fz
831 +          
832 +       endif
833 +    
834         if (nmflag) then
835  
836            drhoidr = drha
# Line 817 | Line 854 | contains
854                 drhojdr*drhojdr*d2frhodrhodrho(atom1)
855   #endif
856         end if
820 –
821 –
857        
858 <      
824 <       if (do_stress) then
825 <          
826 < #ifdef IS_MPI
827 <          id1 = tagRow(atom1)
828 <          id2 = tagColumn(atom2)
829 < #else
830 <          id1 = atom1
831 <          id2 = atom2
832 < #endif
833 <          
834 <          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
835 <
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
839 <             tau_Temp(4) = tau_Temp(4) - d(2) * fx
840 <             tau_Temp(5) = tau_Temp(5) - d(2) * fy
841 <             tau_Temp(6) = tau_Temp(6) - d(2) * fz
842 <             tau_Temp(7) = tau_Temp(7) - d(3) * fx
843 <             tau_Temp(8) = tau_Temp(8) - d(3) * fy
844 <             tau_Temp(9) = tau_Temp(9) - d(3) * fz
845 <
846 <             virial_Temp = virial_Temp + &
847 <                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
848 <
849 <          endif
850 <       endif  
851 <    endif
852 <
853 <    
858 >    endif    
859    end subroutine do_eam_pair
860  
861  

Diff Legend

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