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 631 by chuckv, Thu Jul 17 19:25:51 2003 UTC vs.
Revision 1150 by gezelter, Fri May 7 21:35:05 2004 UTC

# Line 4 | Line 4 | module eam
4    use force_globals
5    use status
6    use atype_module
7 < #ifdef MPI
7 >  use Vector_class
8 > #ifdef IS_MPI
9    use mpiSimulation
10   #endif
11    implicit none
# 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 +  logical, save :: haveRcut = .false.
19  
20 +  character(len = statusMsgSize) :: errMesg
21 +  integer :: eam_err
22 +
23    character(len = 200) :: errMsg
24    character(len=*), parameter :: RoutineName =  "EAM MODULE"
25 + !! Logical that determines if eam arrays should be zeroed
26    logical :: cleanme = .true.
27 +  logical :: nmflag  = .false.
28  
22 –
29  
30    type, private :: EAMtype
31       integer           :: eam_atype      
# Line 49 | Line 55 | module eam
55    real( kind = dp), dimension(:), allocatable :: rho
56  
57    real( kind = dp), dimension(:), allocatable :: dfrhodrho
58 < !  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
58 >  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
59  
60  
61   !! Arrays for MPI storage
62 < #ifdef 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 <
62 > #ifdef IS_MPI
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 71 | 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
89 < !  public :: do_EAM_pair
88 >  public :: setCutoffEAM
89 >  public :: do_eam_pair
90    public :: newEAMtype
91 +  public :: calc_eam_prepair_rho
92 +  public :: calc_eam_preforce_Frho
93 +  public :: clean_EAM
94  
84 –  
85 –
95   contains
96  
97  
# Line 106 | Line 115 | contains
115      integer                                :: alloc_stat
116      integer                                :: current
117      integer,pointer                        :: Matchlist(:) => 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 123 | Line 134 | contains
134      current = EAMList%currentAddition
135      
136  
137 <    !call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),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
142  
132 –
143      ! this is a possible bug, we assume a correspondence between the vector atypes and
144      ! EAMAtypes
145        
# Line 138 | 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 154 | 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 +
177         do i = 1, EAMList%currentAddition
158 –  
159 –          EAMList%EAMParams(i)%eam_rvals(1:EAMList%EAMParams(i)%eam_nr) = &
160 –               real(EAMList%EAMParams(i)%eam_nr,kind=dp)* &
161 –               EAMList%EAMParams(i)%eam_dr
178  
179 <          EAMList%EAMParams(i)%eam_rhovals(1:EAMList%EAMParams(i)%eam_nrho) = &
164 <               real(EAMList%EAMParams(i)%eam_nrho,kind=dp)* &
165 <               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 196 | Line 222 | contains
222                 EAMList%EAMParams(i)%eam_phi_r_pp, &
223                 0.0E0_DP, 0.0E0_DP, 'N')
224         enddo
199 –      
200 –       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
201 –       !! find the smallest rcut
202 –       do i = 2, EAMList%currentAddition
203 –          current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
204 –       end do
225  
226 <       EAM_rcut = current_rcut_max
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 =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
234   !       do i = 1, EAMList%currentAddition
235   !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
236   !       end do
210 –
211 –
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 223 | Line 250 | contains
250    subroutine allocateEAM(status)
251      integer, intent(out) :: status
252  
253 <    integer :: nlocal
227 < #ifdef MPI
253 > #ifdef IS_MPI
254      integer :: nrow
255      integer :: ncol
256   #endif
257      integer :: alloc_stat
258  
259  
260 <    nlocal = getNlocal()
261 <
236 < #ifdef MPI
260 >    status = 0
261 > #ifdef IS_MPI
262      nrow = getNrow(plan_row)
263      ncol = getNcol(plan_col)
264   #endif
# Line 250 | Line 275 | contains
275         status = -1
276         return
277      end if
278 +
279      if (allocated(dfrhodrho)) deallocate(dfrhodrho)
280      allocate(dfrhodrho(nlocal),stat=alloc_stat)
281      if (alloc_stat /= 0) then
282         status = -1
283         return
284      end if
285 +
286 +    if (allocated(d2frhodrhodrho)) deallocate(d2frhodrhodrho)
287 +    allocate(d2frhodrhodrho(nlocal),stat=alloc_stat)
288 +    if (alloc_stat /= 0) then
289 +       status = -1
290 +       return
291 +    end if
292      
293 < #ifdef MPI
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 277 | Line 318 | contains
318         status = -1
319         return
320      end if
321 +    if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
322 +    allocate(d2frhodrhodrho_row(nrow),stat=alloc_stat)
323 +    if (alloc_stat /= 0) then
324 +       status = -1
325 +       return
326 +    end if
327  
328 +
329   ! Now do column arrays
330  
331      if (allocated(frho_col)) deallocate(frho_col)
# Line 298 | Line 346 | contains
346         status = -1
347         return
348      end if
349 <
349 >    if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
350 >    allocate(d2frhodrhodrho_col(ncol),stat=alloc_stat)
351 >    if (alloc_stat /= 0) then
352 >       status = -1
353 >       return
354 >    end if
355 >  
356   #endif
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-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
380   ! clean MPI if needed
381 < #ifdef MPI
381 > #ifdef IS_MPI
382      frho_row = 0.0_dp
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 326 | 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 421 | Line 489 | contains
489      real( kind = dp) :: drho,d2rho
490      integer :: eam_err
491    
492 <    if (cleanme) call clean_EAM
493 <    cleanme = .false.
492 >    integer :: myid_atom1
493 >    integer :: myid_atom2
494  
495 <    call  calc_eam_rho(r,rho_i_at_j,drho,d2rho,atom1)
495 > ! check to see if we need to be cleaned at the start of a force loop
496 >      
497  
498 < #ifdef  MPI
499 <    rho_col(atom2) = rho_col(atom2) + rho_i_at_j
498 >
499 >
500 > #ifdef IS_MPI
501 >    myid_atom1 = atid_Row(atom1)
502 >    myid_atom2 = atid_Col(atom2)
503   #else
504 <    rho(atom2) = rho(atom2) + rho_i_at_j
504 >    myid_atom1 = atid(atom1)
505 >    myid_atom2 = atid(atom2)
506   #endif
507  
508 <    call calc_eam_rho(r,rho_j_at_i,drho,d2rho,atom2)
508 >    if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
509  
510 < #ifdef  MPI
511 <    rho_row(atom1) = rho_row(atom1) + rho_j_at_i
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, &
515 >            EAMList%EAMParams(myid_atom1)%eam_rho_r_pp, &
516 >            r, rho_i_at_j,drho,d2rho)
517 >
518 >
519 >      
520 > #ifdef  IS_MPI
521 >       rho_col(atom2) = rho_col(atom2) + rho_i_at_j
522   #else
523 <    rho(atom1) = rho(atom1) + rho_j_at_i
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
529 +          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
530 +               EAMList%EAMParams(myid_atom2)%eam_rvals, &
531 +               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
532 +               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
533 +               r, rho_j_at_i,drho,d2rho)
534 +
535 +    
536 +      
537 +      
538 + #ifdef  IS_MPI
539 +          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
540 + #else
541 +          rho(atom1) = rho(atom1) + rho_j_at_i
542 + #endif
543 +       endif
544 +
545 +
546 +
547 +
548 +
549 +
550    end subroutine calc_eam_prepair_rho
551  
552 <  !! Calculate the functional F(rho) for all atoms
553 <  subroutine calc_eam_prepair_Frho(nlocal,pot)
552 >
553 >
554 >
555 >  !! Calculate the functional F(rho) for all local atoms
556 >  subroutine calc_eam_preforce_Frho(nlocal,pot)
557      integer :: nlocal
558      real(kind=dp) :: pot
559      integer :: i,j
560 +    integer :: atom
561      real(kind=dp) :: U,U1,U2
562      integer :: atype1
563 <    ! reset clean forces to be true at top of calc rho.
564 <    cleanme = .true.
563 >    integer :: me
564 >    integer :: n_rho_points
565  
566 < !! Scatter the electron density in pre-pair
567 < #ifdef MPI
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  
468 –   do i = 1, nlocal
469 –      call calc_eam_frho(rho(i),u,u1,u2,atype1)
470 –      frho(i) = u
471 –      dfrhodrho(i) = u1
472 – !      d2frhodrhodrho(i) = u2
473 –      pot = pot + u
474 –   enddo
584  
585 < #ifdef MPI
586 <    !! communicate f(rho) and derivatives
585 >
586 > !! Calculate F(rho) and derivative
587 >    do atom = 1, nlocal
588 >       me = atid(atom)
589 >       n_rho_points = EAMList%EAMParams(me)%eam_nrho
590 >       !  Check to see that the density is not greater than the larges rho we have calculated
591 >       if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
592 >          call eam_splint(n_rho_points, &
593 >               EAMList%EAMParams(me)%eam_rhovals, &
594 >               EAMList%EAMParams(me)%eam_f_rho, &
595 >               EAMList%EAMParams(me)%eam_f_rho_pp, &
596 >               rho(atom), & ! Actual Rho
597 >               u, u1, u2)
598 >       else
599 >          ! Calculate F(rho with the largest available rho value
600 >          call eam_splint(n_rho_points, &
601 >               EAMList%EAMParams(me)%eam_rhovals, &
602 >               EAMList%EAMParams(me)%eam_f_rho, &
603 >               EAMList%EAMParams(me)%eam_f_rho_pp, &
604 >               EAMList%EAMParams(me)%eam_rhovals(n_rho_points), & ! Largest rho
605 >               u,u1,u2)
606 >       end if
607 >
608 >
609 >       frho(atom) = u
610 >       dfrhodrho(atom) = u1
611 >       d2frhodrhodrho(atom) = u2
612 >       pot = pot + u
613 >
614 >    enddo
615 >
616 >  
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)
621      if (eam_err /=  0) then
622         call handleError("cal_eam_forces()","MPI gather frho_row failure")
# Line 491 | Line 633 | contains
633      if (eam_err /=  0) then
634         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
635      endif
636 +
637 +
638  
639 +
640 +
641      if (nmflag) then
642         call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_row)
643         call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_col)
644      endif
645   #endif
646  
647 +  
648 +  end subroutine calc_eam_preforce_Frho
649  
502 –  end subroutine calc_eam_prepair_Frho
650  
651  
652  
653 <
654 <  subroutine calc_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress)
655 < !Arguments    
653 >  !! Does EAM pairwise Force calculation.  
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 <
663 >    
664      real( kind = dp ) :: drdx,drdy,drdz
665 +    real( kind = dp ) :: d2
666      real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
667      real( kind = dp ) :: rha,drha,d2rha, dpha
668      real( kind = dp ) :: rhb,drhb,d2rhb, dphb
# Line 524 | Line 673 | contains
673      real( kind = dp ) :: d2rhojdrdr
674      real( kind = dp ) :: Fx,Fy,Fz
675      real( kind = dp ) :: r,d2pha,phb,d2phb
676 +
677      integer :: id1,id2
678 <    integer  :: atype1,atype2
678 >    integer  :: mytype_atom1
679 >    integer  :: mytype_atom2
680  
530 –
681   !Local Variables
682      
683 +   ! write(*,*) "Frho: ", Frho(atom1)
684  
534 –    
685      phab = 0.0E0_DP
686      dvpdr = 0.0E0_DP
687      d2vpdrdr = 0.0E0_DP
538 –      
688  
689      if (rij .lt. EAM_rcut) then
690   #ifdef IS_MPI
691   !!!!! FIX ME
692 <       atype1 = atid_row(atom1)
692 >       mytype_atom1 = atid_row(atom1)
693   #else
694 <       atype1 = atid(atom1)
694 >       mytype_atom1 = atid(atom1)
695   #endif
696 <
696 >      
697         drdx = d(1)/rij
698         drdy = d(2)/rij
699         drdz = d(3)/rij
700        
701 <      
702 <       call calc_eam_rho(r, rha, drha, d2rha, atype1)
703 <       call calc_eam_phi(r, pha, dpha, d2pha, atype1)
704 <  !     rci = eam_rcut(eam_atype_map(atom1))
701 >
702 >       call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
703 >               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
704 >               EAMList%EAMParams(mytype_atom1)%eam_rho_r, &
705 >               EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, &
706 >               rij, rha,drha,d2rha)
707 >
708 >       !! Calculate Phi(r) for atom1.
709 >       call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
710 >               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
711 >               EAMList%EAMParams(mytype_atom1)%eam_phi_r, &
712 >               EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, &
713 >               rij, pha,dpha,d2pha)
714 >
715 >
716 > ! get cutoff for atom 1
717 >       rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
718   #ifdef IS_MPI
719 <       atype2 = atid_col(atom2)
719 >       mytype_atom2 = atid_col(atom2)
720   #else
721 <       atype2 = atid(atom2)
721 >       mytype_atom2 = atid(atom2)
722   #endif
561 –      
562 –       call calc_eam_rho(r, rhb, drhb, d2rhb, atype2)
563 –       call calc_eam_phi(r, phb, dphb, d2phb, atype2)
564 –  !     rcj = eam_rcut(eam_atype_map(atype2))
723  
724 <       if (r.lt.rci) then
724 >       ! Calculate rho,drho and d2rho for atom1
725 >       call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
726 >               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
727 >               EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
728 >               EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
729 >               rij, rhb,drhb,d2rhb)
730 >
731 >       !! Calculate Phi(r) for atom2.
732 >       call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
733 >               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
734 >               EAMList%EAMParams(mytype_atom2)%eam_phi_r, &
735 >               EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, &
736 >               rij, phb,dphb,d2phb)
737 >
738 >
739 > ! get type specific cutoff for atom 2
740 >       rcj = EAMList%EAMParams(mytype_atom1)%eam_rcut
741 >
742 >
743 >
744 >       if (rij.lt.rci) then
745            phab = phab + 0.5E0_DP*(rhb/rha)*pha
746            dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
747                 pha*((drhb/rha) - (rhb*drha/rha/rha)))
# Line 574 | Line 752 | contains
752         endif
753        
754  
755 <       if (r.lt.rcj) then
755 >       if (rij.lt.rcj) then
756            phab = phab + 0.5E0_DP*(rha/rhb)*phb
757            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
758                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)))
# Line 591 | Line 769 | contains
769         d2rhojdrdr = d2rhb
770  
771  
772 < #ifdef MPI
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 606 | Line 784 | contains
784         fz = dudr * drdz
785  
786  
787 < #ifdef MPI
787 > #ifdef IS_MPI
788         if (do_pot) then
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
796         f_Row(3,atom1) = f_Row(3,atom1) + fz
797 <
797 >      
798         f_Col(1,atom2) = f_Col(1,atom2) - fx
799         f_Col(2,atom2) = f_Col(2,atom2) - fy
800         f_Col(3,atom2) = f_Col(3,atom2) - fz
801   #else
623 –       if(do_pot) pot = pot + phab
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
811 <
811 >      
812         f(1,atom2) = f(1,atom2) - fx
813         f(2,atom2) = f(2,atom2) - fy
814         f(3,atom2) = f(3,atom2) - fz
815   #endif
816 +      
817 +       if (nmflag) then
818  
819 +          drhoidr = drha
820 +          drhojdr = drhb
821 +          d2rhoidrdr = d2rha
822 +          d2rhojdrdr = d2rhb
823  
824 + #ifdef IS_MPI
825 +          d2 = d2vpdrdr + &
826 +               d2rhoidrdr*dfrhodrho_col(atom2) + &
827 +               d2rhojdrdr*dfrhodrho_row(atom1) + &
828 +               drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
829 +               drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
830 +              
831 + #else
832  
833 +          d2 = d2vpdrdr + &
834 +               d2rhoidrdr*dfrhodrho(atom2) + &
835 +               d2rhojdrdr*dfrhodrho(atom1) + &
836 +               drhoidr*drhoidr*d2frhodrhodrho(atom2) + &
837 +               drhojdr*drhojdr*d2frhodrhodrho(atom1)
838 + #endif
839 +       end if
840  
637 –       if (do_stress) then
841  
842 < #ifdef MPI
842 >      
843 >      
844 >       if (do_stress) then
845 >          
846 > #ifdef IS_MPI
847            id1 = tagRow(atom1)
848            id2 = tagColumn(atom2)
849   #else
850            id1 = atom1
851            id2 = atom2
852   #endif
853 <
853 >          
854            if (molMembershipList(id1) .ne. molMembershipList(id2)) then
855 <            
856 <             tau_Temp(1) = tau_Temp(1) + fx * d(1)
857 <             tau_Temp(2) = tau_Temp(2) + fx * d(2)
858 <             tau_Temp(3) = tau_Temp(3) + fx * d(3)
859 <             tau_Temp(4) = tau_Temp(4) + fy * d(1)
860 <             tau_Temp(5) = tau_Temp(5) + fy * d(2)
861 <             tau_Temp(6) = tau_Temp(6) + fy * d(3)
862 <             tau_Temp(7) = tau_Temp(7) + fz * d(1)
863 <             tau_Temp(8) = tau_Temp(8) + fz * d(2)
864 <             tau_Temp(9) = tau_Temp(9) + fz * d(3)
855 >
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
859 >             tau_Temp(4) = tau_Temp(4) - d(2) * fx
860 >             tau_Temp(5) = tau_Temp(5) - d(2) * fy
861 >             tau_Temp(6) = tau_Temp(6) - d(2) * fz
862 >             tau_Temp(7) = tau_Temp(7) - d(3) * fx
863 >             tau_Temp(8) = tau_Temp(8) - d(3) * fy
864 >             tau_Temp(9) = tau_Temp(9) - d(3) * fz
865 >
866               virial_Temp = virial_Temp + &
867                    (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
868  
869            endif
870 <       endif
663 <      
870 >       endif  
871      endif
872  
873 +    
874 +  end subroutine do_eam_pair
875  
667 – end subroutine calc_eam_pair
876  
669 – !!$subroutine calc_eam_rho(r, rho, drho, d2rho, atype)
670 – !!$
671 – !!$  !  include 'headers/sizes.h'
672 – !!$
673 – !!$
674 – !!$integer atype, etype, number_r
675 – !!$real( kind = DP )  :: r, rho, drho, d2rho
676 – !!$integer :: i
677 – !!$
678 – !!$
679 – !!$etype = eam_atype_map(atype)
680 – !!$
681 – !!$if (r.lt.eam_rcut(etype)) then
682 – !!$number_r = eam_nr(etype)
683 – !!$call eam_splint(etype, number_r, eam_rvals, eam_rho_r, &
684 – !!$   eam_rho_r_pp, r, rho, drho, d2rho)
685 – !!$else
686 – !!$rho = 0.0E0_DP
687 – !!$drho = 0.0E0_DP
688 – !!$d2rho = 0.0E0_DP
689 – !!$endif
690 – !!$
691 – !!$return
692 – !!$end subroutine calc_eam_rho
693 – !!$
694 – !!$subroutine calc_eam_frho(dens, u, u1, u2, atype)
695 – !!$
696 – !!$  ! include 'headers/sizes.h'
697 – !!$
698 – !!$integer atype, etype, number_rho
699 – !!$real( kind = DP ) :: dens, u, u1, u2
700 – !!$real( kind = DP ) :: rho_vals
701 – !!$
702 – !!$etype = eam_atype_map(atype)
703 – !!$number_rho = eam_nrho(etype)
704 – !!$if (dens.lt.eam_rhovals(number_rho, etype)) then
705 – !!$call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
706 – !!$   eam_f_rho_pp, dens, u, u1, u2)
707 – !!$else
708 – !!$rho_vals = eam_rhovals(number_rho,etype)
709 – !!$call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
710 – !!$   eam_f_rho_pp, rho_vals, u, u1, u2)
711 – !!$endif
712 – !!$
713 – !!$return
714 – !!$end subroutine calc_eam_frho
715 – !!$
716 – !!$subroutine calc_eam_phi(r, phi, dphi, d2phi, atype)
717 – !!$
718 – !!$
719 – !!$
720 – !!$
721 – !!$integer atype, etype, number_r
722 – !!$real( kind = DP ) :: r, phi, dphi, d2phi
723 – !!$
724 – !!$etype = eam_atype_map(atype)
725 – !!$
726 – !!$if (r.lt.eam_rcut(etype)) then
727 – !!$number_r = eam_nr(etype)
728 – !!$call eam_splint(etype, number_r, eam_rvals, eam_phi_r, &
729 – !!$   eam_phi_r_pp, r, phi, dphi, d2phi)
730 – !!$else
731 – !!$phi = 0.0E0_DP
732 – !!$dphi = 0.0E0_DP
733 – !!$d2phi = 0.0E0_DP
734 – !!$endif
735 – !!$
736 – !!$return
737 – !!$end subroutine calc_eam_phi
738 –
739 –
877    subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
878  
742 –
879      integer :: atype, nx, j
880      real( kind = DP ), dimension(:) :: xa
881      real( kind = DP ), dimension(:) :: ya
882      real( kind = DP ), dimension(:) :: yppa
883 <    real( kind = DP ) :: x, y, dy, d2y
883 >    real( kind = DP ) :: x, y
884 >    real( kind = DP ) :: dy, d2y
885      real( kind = DP ) :: del, h, a, b, c, d
886 +    integer :: pp_arraySize
887  
888 <
751 <
752 <
753 <    
888 >
889      ! this spline code assumes that the x points are equally spaced
890      ! do not attempt to use this code if they are not.
891      
892      
893      ! find the closest point with a value below our own:
894 <    j = FLOOR(dble(nx-1) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
894 >    j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
895  
896      ! check to make sure we're inside the spline range:
897      if ((j.gt.nx).or.(j.lt.1)) then
898 <       write(default_error,*) "EAM_splint: x is outside bounds of spline"
898 >       write(errMSG,*) "EAM_splint: x is outside bounds of spline"
899 >       call handleError(routineName,errMSG)
900      endif
901      ! check to make sure we haven't screwed up the calculation of j:
902      if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
903         if (j.ne.nx) then
904 <        write(default_error,*) "EAM_splint: x is outside bounding range"
904 >        write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
905 >       call handleError(routineName,errMSG)
906         endif
907      endif
908  
# Line 778 | Line 915 | end subroutine calc_eam_pair
915      d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
916      
917      y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
918 <    
919 <    dy = (ya(j+1)-ya(j))/h &
920 <         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
921 <         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
922 <    
923 <    d2y = a*yppa(j) + b*yppa(j+1)
918 >  
919 >       dy = (ya(j+1)-ya(j))/h &
920 >            - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
921 >            + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
922 >  
923 >  
924 >       d2y = a*yppa(j) + b*yppa(j+1)
925 >  
926  
927    end subroutine eam_splint
928  
929 +
930    subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
931  
792 –  
932  
794 –
933      ! yp1 and ypn are the first derivatives of y at the two endpoints
934      ! if boundary is 'L' the lower derivative is used
935      ! if boundary is 'U' the upper derivative is used
# Line 805 | Line 943 | end subroutine calc_eam_pair
943      real( kind = DP ), dimension(:)        :: yppa
944      real( kind = DP ), dimension(size(xa)) :: u
945      real( kind = DP ) :: yp1,ypn,un,qn,sig,p
946 <    character boundary
946 >    character(len=*) :: boundary
947      
948 <    
948 >    ! make sure the sizes match
949 >    if ((nx /= size(xa)) .or. (nx /= size(ya))) then
950 >       call handleWarning("EAM_SPLINE","Array size mismatch")
951 >    end if
952 >
953      if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
954           (boundary.eq.'b').or.(boundary.eq.'B')) then
955         yppa(1) = -0.5E0_DP

Diff Legend

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