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 669 by chuckv, Thu Aug 7 00:47:33 2003 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 45 | Line 51 | module eam
51  
52  
53    !! Arrays for derivatives used in force calculation
54 <  real( kind = dp), dimension(:), allocatable :: frho
55 <  real( kind = dp), dimension(:), allocatable :: rho
54 >  real( kind = dp),save, dimension(:), allocatable :: frho
55 >  real( kind = dp),save, dimension(:), allocatable :: rho
56  
57 <  real( kind = dp), dimension(:), allocatable :: dfrhodrho
58 < !  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
57 >  real( kind = dp),save, dimension(:), allocatable :: dfrhodrho
58 >  real( kind = dp),save, 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 :: d2frhodrhodrho_col
70 >  real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
71   #endif
72  
73    type, private :: EAMTypeList
# Line 77 | Line 84 | module eam
84  
85  
86    public :: init_EAM_FF
87 < !  public :: EAM_new_rcut
88 < !  public :: do_EAM_pair
87 >  public :: setCutoffEAM
88 >  public :: do_eam_pair
89    public :: newEAMtype
90 <
90 >  public :: calc_eam_prepair_rho
91 >  public :: calc_eam_preforce_Frho
92    
93  
94   contains
# Line 106 | Line 114 | contains
114      integer                                :: alloc_stat
115      integer                                :: current
116      integer,pointer                        :: Matchlist(:) => 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 123 | Line 133 | contains
133      current = EAMList%currentAddition
134      
135  
136 <    !call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),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
141  
132 –
142      ! this is a possible bug, we assume a correspondence between the vector atypes and
143      ! EAMAtypes
144        
# Line 138 | 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 154 | Line 164 | contains
164      integer :: alloc_stat
165      integer :: number_r, number_rho
166  
167 +    if (EAMList%currentAddition == 0) then
168 +       call handleError("init_EAM_FF","No members in EAMList")
169 +       status = -1
170 +       return
171 +    end if
172 +
173 +
174 +
175         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
176  
177 <          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
177 > ! Build array of r values
178  
179 +          do j = 1,EAMList%EAMParams(i)%eam_nr
180 +             EAMList%EAMParams(i)%eam_rvals(j) = &
181 +                  real(j-1,kind=dp)* &
182 +                  EAMList%EAMParams(i)%eam_dr
183 +              end do
184 + ! Build array of rho values
185 +          do j = 1,EAMList%EAMParams(i)%eam_nrho
186 +             EAMList%EAMParams(i)%eam_rhovals(j) = &
187 +                  real(j-1,kind=dp)* &
188 +                  EAMList%EAMParams(i)%eam_drho
189 +          end do
190            ! convert from eV to kcal / mol:
191            EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
192  
193            ! precompute the pair potential and get it into kcal / mol:
194 <          EAMList%EAMParams(i)%eam_phi_r = 0.0E0_DP
194 >          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
195            do j = 2, EAMList%EAMParams(i)%eam_nr
196               EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
197               EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
198            enddo
199         end do
200 +      
201  
202         do i = 1,  EAMList%currentAddition
203            number_r   = EAMList%EAMParams(i)%eam_nr
# Line 196 | Line 220 | contains
220                 EAMList%EAMParams(i)%eam_phi_r_pp, &
221                 0.0E0_DP, 0.0E0_DP, 'N')
222         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
223  
224 <       EAM_rcut = current_rcut_max
224 >
225 > !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
226 >       !! find the smallest rcut for any eam atype
227 > !       do i = 2, EAMList%currentAddition
228 > !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
229 > !       end do
230 >
231 > !       EAM_rcut = current_rcut_max
232 > !       EAM_rcut_orig = current_rcut_max
233   !       do i = 1, EAMList%currentAddition
234   !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
235   !       end do
# Line 224 | Line 250 | contains
250      integer, intent(out) :: status
251  
252      integer :: nlocal
253 < #ifdef MPI
253 > #ifdef IS_MPI
254      integer :: nrow
255      integer :: ncol
256   #endif
# Line 233 | Line 259 | contains
259  
260      nlocal = getNlocal()
261  
262 < #ifdef MPI
262 > #ifdef IS_MPI
263      nrow = getNrow(plan_row)
264      ncol = getNcol(plan_col)
265   #endif
# Line 250 | Line 276 | contains
276         status = -1
277         return
278      end if
279 +
280      if (allocated(dfrhodrho)) deallocate(dfrhodrho)
281      allocate(dfrhodrho(nlocal),stat=alloc_stat)
282      if (alloc_stat /= 0) then
283         status = -1
284         return
285      end if
286 +
287 +    if (allocated(d2frhodrhodrho)) deallocate(d2frhodrhodrho)
288 +    allocate(d2frhodrhodrho(nlocal),stat=alloc_stat)
289 +    if (alloc_stat /= 0) then
290 +       status = -1
291 +       return
292 +    end if
293      
294 < #ifdef MPI
294 > #ifdef IS_MPI
295  
296      if (allocated(frho_row)) deallocate(frho_row)
297      allocate(frho_row(nrow),stat=alloc_stat)
# Line 277 | Line 311 | contains
311         status = -1
312         return
313      end if
314 +    if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
315 +    allocate(d2frhodrhodrho_row(nrow),stat=alloc_stat)
316 +    if (alloc_stat /= 0) then
317 +       status = -1
318 +       return
319 +    end if
320  
321 +
322   ! Now do column arrays
323  
324      if (allocated(frho_col)) deallocate(frho_col)
# Line 298 | Line 339 | contains
339         status = -1
340         return
341      end if
342 <
342 >    if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
343 >    allocate(d2frhodrhodrho_col(ncol),stat=alloc_stat)
344 >    if (alloc_stat /= 0) then
345 >       status = -1
346 >       return
347 >    end if
348 >  
349   #endif
350  
351    end subroutine allocateEAM
352  
353 + !! C sets rcut to be the largest cutoff of any atype
354 + !! present in this simulation. Doesn't include all atypes
355 + !! sim knows about, just those in the simulation.
356 +  subroutine setCutoffEAM(rcut, status)
357 +    real(kind=dp) :: rcut
358 +    integer :: status
359 +    status = 0
360  
361 <  subroutine clean_EAM()
361 >    EAM_rcut = rcut
362  
363 < ! clean non-MPI first
363 >  end subroutine setCutoffEAM
364 >
365 >
366 >
367 >  subroutine clean_EAM()
368 >   ! clean non-IS_MPI first
369      frho = 0.0_dp
370      rho  = 0.0_dp
371      dfrhodrho = 0.0_dp
372   ! clean MPI if needed
373 < #ifdef MPI
373 > #ifdef IS_MPI
374      frho_row = 0.0_dp
375      frho_col = 0.0_dp
376      rho_row  = 0.0_dp
# Line 326 | Line 385 | contains
385    subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
386      integer, intent(in)          :: eam_n_rho
387      integer, intent(in)          :: eam_n_r
388 <    type (EAMType), pointer      :: thisEAMType
388 >    type (EAMType)               :: thisEAMType
389      integer, optional   :: stat
390      integer             :: alloc_stat
391  
# Line 421 | Line 480 | contains
480      real( kind = dp) :: drho,d2rho
481      integer :: eam_err
482    
483 <    if (cleanme) call clean_EAM
484 <    cleanme = .false.
483 >    integer :: myid_atom1
484 >    integer :: myid_atom2
485  
486 <    call  calc_eam_rho(r,rho_i_at_j,drho,d2rho,atom1)
486 > ! 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
492 >      
493  
494 < #ifdef  MPI
495 <    rho_col(atom2) = rho_col(atom2) + rho_i_at_j
494 >    
495 >
496 > #ifdef IS_MPI
497 >    myid_atom1 = atid_Row(atom1)
498 >    myid_atom2 = atid_Col(atom2)
499   #else
500 <    rho(atom2) = rho(atom2) + rho_i_at_j
500 >    myid_atom1 = atid(atom1)
501 >    myid_atom2 = atid(atom2)
502   #endif
503  
504 <    call calc_eam_rho(r,rho_j_at_i,drho,d2rho,atom2)
504 >    if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
505 >
506  
507 < #ifdef  MPI
508 <    rho_row(atom1) = rho_row(atom1) + rho_j_at_i
507 >
508 >       call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
509 >            EAMList%EAMParams(myid_atom1)%eam_rvals, &
510 >            EAMList%EAMParams(myid_atom1)%eam_rho_r, &
511 >            EAMList%EAMParams(myid_atom1)%eam_rho_r_pp, &
512 >            r, rho_i_at_j,drho,d2rho)
513 >
514 >
515 >      
516 > #ifdef  IS_MPI
517 >       rho_col(atom2) = rho_col(atom2) + rho_i_at_j
518   #else
519 <    rho(atom1) = rho(atom1) + rho_j_at_i
519 >       rho(atom2) = rho(atom2) + rho_i_at_j
520   #endif
521 +       endif
522 +
523 +       if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
524 +          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
525 +               EAMList%EAMParams(myid_atom2)%eam_rvals, &
526 +               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
527 +               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
528 +               r, rho_j_at_i,drho,d2rho)
529 +
530 +
531 +      
532 +      
533 + #ifdef  IS_MPI
534 +          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
535 + #else
536 +          rho(atom1) = rho(atom1) + rho_j_at_i
537 + #endif
538 +       endif
539 +
540 +    
541 +
542    end subroutine calc_eam_prepair_rho
543  
544 <  !! Calculate the functional F(rho) for all atoms
545 <  subroutine calc_eam_prepair_Frho(nlocal,pot)
544 >
545 >
546 >
547 >  !! Calculate the functional F(rho) for all local atoms
548 >  subroutine calc_eam_preforce_Frho(nlocal,pot)
549      integer :: nlocal
550      real(kind=dp) :: pot
551      integer :: i,j
552 +    integer :: atom
553      real(kind=dp) :: U,U1,U2
554      integer :: atype1
555 +    integer :: me
556 +    integer :: n_rho_points
557      ! reset clean forces to be true at top of calc rho.
558      cleanme = .true.
559 <
560 < !! Scatter the electron density in pre-pair
561 < #ifdef MPI
559 >  
560 > !! Scatter the electron density from  pre-pair calculation back to local atoms
561 > #ifdef IS_MPI
562      call scatter(rho_row,rho,plan_row,eam_err)
563      if (eam_err /= 0 ) then
564        write(errMsg,*) " Error scattering rho_row into rho"
# Line 465 | Line 571 | contains
571     endif
572   #endif
573  
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
574  
575 < #ifdef MPI
576 <    !! communicate f(rho) and derivatives
575 > !! Calculate F(rho) and derivative
576 >    do atom = 1, nlocal
577 >       me = atid(atom)
578 >       n_rho_points = EAMList%EAMParams(me)%eam_nrho
579 >       !  Check to see that the density is not greater than the larges rho we have calculated
580 >       if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
581 >          call eam_splint(n_rho_points, &
582 >               EAMList%EAMParams(me)%eam_rhovals, &
583 >               EAMList%EAMParams(me)%eam_f_rho, &
584 >               EAMList%EAMParams(me)%eam_f_rho_pp, &
585 >               rho(atom), & ! Actual Rho
586 >               u, u1, u2)
587 >       else
588 >          ! Calculate F(rho with the largest available rho value
589 >          call eam_splint(n_rho_points, &
590 >               EAMList%EAMParams(me)%eam_rhovals, &
591 >               EAMList%EAMParams(me)%eam_f_rho, &
592 >               EAMList%EAMParams(me)%eam_f_rho_pp, &
593 >               EAMList%EAMParams(me)%eam_rhovals(n_rho_points), & ! Largest rho
594 >               u,u1,u2)
595 >       end if
596 >
597 >
598 >       frho(i) = u
599 >       dfrhodrho(i) = u1
600 >       d2frhodrhodrho(i) = u2
601 >       pot = pot + u
602 >    enddo
603 >
604 >  
605 >
606 > #ifdef IS_MPI
607 >    !! communicate f(rho) and derivatives back into row and column arrays
608      call gather(frho,frho_row,plan_row, eam_err)
609      if (eam_err /=  0) then
610         call handleError("cal_eam_forces()","MPI gather frho_row failure")
# Line 492 | Line 622 | contains
622         call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
623      endif
624  
625 +
626 +
627 +
628 +
629      if (nmflag) then
630         call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_row)
631         call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_col)
632      endif
633   #endif
634  
635 +  end subroutine calc_eam_preforce_Frho
636  
502 –  end subroutine calc_eam_prepair_Frho
637  
638  
639  
640 <
641 <  subroutine calc_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress)
642 < !Arguments    
640 >  !! Does EAM pairwise Force calculation.  
641 >  subroutine do_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress)
642 >    !Arguments    
643      integer, intent(in) ::  atom1, atom2
644      real( kind = dp ), intent(in) :: rij, r2
645      real( kind = dp ) :: pot
646      real( kind = dp ), dimension(3,getNlocal()) :: f
647      real( kind = dp ), intent(in), dimension(3) :: d
648      logical, intent(in) :: do_pot, do_stress
649 <
649 >    
650      real( kind = dp ) :: drdx,drdy,drdz
651 +    real( kind = dp ) :: d2
652      real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
653      real( kind = dp ) :: rha,drha,d2rha, dpha
654      real( kind = dp ) :: rhb,drhb,d2rhb, dphb
# Line 524 | Line 659 | contains
659      real( kind = dp ) :: d2rhojdrdr
660      real( kind = dp ) :: Fx,Fy,Fz
661      real( kind = dp ) :: r,d2pha,phb,d2phb
662 +
663      integer :: id1,id2
664 <    integer  :: atype1,atype2
664 >    integer  :: mytype_atom1
665 >    integer  :: mytype_atom2
666  
667  
668   !Local Variables
669      
670  
671 <    
671 >
672      phab = 0.0E0_DP
673      dvpdr = 0.0E0_DP
674      d2vpdrdr = 0.0E0_DP
538 –      
675  
676      if (rij .lt. EAM_rcut) then
677   #ifdef IS_MPI
678   !!!!! FIX ME
679 <       atype1 = atid_row(atom1)
679 >       mytype_atom1 = atid_row(atom1)
680   #else
681 <       atype1 = atid(atom1)
681 >       mytype_atom1 = atid(atom1)
682   #endif
683 <
683 >      
684         drdx = d(1)/rij
685         drdy = d(2)/rij
686         drdz = d(3)/rij
687        
688 <      
689 <       call calc_eam_rho(r, rha, drha, d2rha, atype1)
690 <       call calc_eam_phi(r, pha, dpha, d2pha, atype1)
691 <  !     rci = eam_rcut(eam_atype_map(atom1))
688 >
689 >       call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
690 >               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
691 >               EAMList%EAMParams(mytype_atom1)%eam_rho_r, &
692 >               EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, &
693 >               rij, rha,drha,d2rha)
694 >
695 >       !! Calculate Phi(r) for atom1.
696 >       call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
697 >               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
698 >               EAMList%EAMParams(mytype_atom1)%eam_phi_r, &
699 >               EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, &
700 >               rij, pha,dpha,d2pha)
701 >
702 >
703 > ! get cutoff for atom 1
704 >       rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
705   #ifdef IS_MPI
706 <       atype2 = atid_col(atom2)
706 >       mytype_atom2 = atid_col(atom2)
707   #else
708 <       atype2 = atid(atom2)
708 >       mytype_atom2 = atid(atom2)
709   #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))
710  
711 <       if (r.lt.rci) then
711 >       ! Calculate rho,drho and d2rho for atom1
712 >       call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
713 >               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
714 >               EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
715 >               EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
716 >               rij, rhb,drhb,d2rhb)
717 >
718 >       !! Calculate Phi(r) for atom2.
719 >       call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
720 >               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
721 >               EAMList%EAMParams(mytype_atom2)%eam_phi_r, &
722 >               EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, &
723 >               rij, phb,dphb,d2phb)
724 >
725 >
726 > ! get type specific cutoff for atom 2
727 >       rcj = EAMList%EAMParams(mytype_atom1)%eam_rcut
728 >
729 >
730 >
731 >       if (rij.lt.rci) then
732            phab = phab + 0.5E0_DP*(rhb/rha)*pha
733            dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
734                 pha*((drhb/rha) - (rhb*drha/rha/rha)))
# Line 574 | Line 739 | contains
739         endif
740        
741  
742 <       if (r.lt.rcj) then
742 >       if (rij.lt.rcj) then
743            phab = phab + 0.5E0_DP*(rha/rhb)*phb
744            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
745                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)))
# Line 591 | Line 756 | contains
756         d2rhojdrdr = d2rhb
757  
758  
759 < #ifdef MPI
759 > #ifdef IS_MPI
760         dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
761              + dvpdr
762 <    
762 >
763   #else
764         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
765              + dvpdr
766 <
766 >      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
767   #endif
768  
769         fx = dudr * drdx
# Line 606 | Line 771 | contains
771         fz = dudr * drdz
772  
773  
774 < #ifdef MPI
774 > #ifdef IS_MPI
775         if (do_pot) then
776            pot_Row(atom1) = pot_Row(atom1) + phab*0.5
777            pot_Col(atom2) = pot_Col(atom2) + phab*0.5
# Line 615 | Line 780 | contains
780         f_Row(1,atom1) = f_Row(1,atom1) + fx
781         f_Row(2,atom1) = f_Row(2,atom1) + fy
782         f_Row(3,atom1) = f_Row(3,atom1) + fz
783 <
783 >      
784         f_Col(1,atom2) = f_Col(1,atom2) - fx
785         f_Col(2,atom2) = f_Col(2,atom2) - fy
786         f_Col(3,atom2) = f_Col(3,atom2) - fz
787   #else
623 –       if(do_pot) pot = pot + phab
788  
789 +       if(do_pot) then
790 +          pot = pot + phab
791 +       end if
792 +
793         f(1,atom1) = f(1,atom1) + fx
794         f(2,atom1) = f(2,atom1) + fy
795         f(3,atom1) = f(3,atom1) + fz
796 <
796 >      
797         f(1,atom2) = f(1,atom2) - fx
798         f(2,atom2) = f(2,atom2) - fy
799         f(3,atom2) = f(3,atom2) - fz
800   #endif
801 +      
802 +       if (nmflag) then
803  
804 +          drhoidr = drha
805 +          drhojdr = drhb
806 +          d2rhoidrdr = d2rha
807 +          d2rhojdrdr = d2rhb
808  
809 + #ifdef IS_MPI
810 +          d2 = d2vpdrdr + &
811 +               d2rhoidrdr*dfrhodrho_col(atom2) + &
812 +               d2rhojdrdr*dfrhodrho_row(atom1) + &
813 +               drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
814 +               drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
815 +              
816 + #else
817  
818 +          d2 = d2vpdrdr + &
819 +               d2rhoidrdr*dfrhodrho(atom2) + &
820 +               d2rhojdrdr*dfrhodrho(atom1) + &
821 +               drhoidr*drhoidr*d2frhodrhodrho(atom2) + &
822 +               drhojdr*drhojdr*d2frhodrhodrho(atom1)
823 + #endif
824 +       end if
825  
637 –       if (do_stress) then
826  
827 < #ifdef MPI
827 >      
828 >      
829 >       if (do_stress) then
830 >          
831 > #ifdef IS_MPI
832            id1 = tagRow(atom1)
833            id2 = tagColumn(atom2)
834   #else
835            id1 = atom1
836            id2 = atom2
837   #endif
838 <
838 >          
839            if (molMembershipList(id1) .ne. molMembershipList(id2)) then
840              
841 <             tau_Temp(1) = tau_Temp(1) + fx * d(1)
842 <             tau_Temp(2) = tau_Temp(2) + fx * d(2)
843 <             tau_Temp(3) = tau_Temp(3) + fx * d(3)
844 <             tau_Temp(4) = tau_Temp(4) + fy * d(1)
845 <             tau_Temp(5) = tau_Temp(5) + fy * d(2)
846 <             tau_Temp(6) = tau_Temp(6) + fy * d(3)
847 <             tau_Temp(7) = tau_Temp(7) + fz * d(1)
848 <             tau_Temp(8) = tau_Temp(8) + fz * d(2)
849 <             tau_Temp(9) = tau_Temp(9) + fz * d(3)
841 >
842 >
843 >
844 >             tau_Temp(1) = tau_Temp(1) - d(1) * fx
845 >             tau_Temp(2) = tau_Temp(2) - d(1) * fy
846 >             tau_Temp(3) = tau_Temp(3) - d(1) * fz
847 >             tau_Temp(4) = tau_Temp(4) - d(2) * fx
848 >             tau_Temp(5) = tau_Temp(5) - d(2) * fy
849 >             tau_Temp(6) = tau_Temp(6) - d(2) * fz
850 >             tau_Temp(7) = tau_Temp(7) - d(3) * fx
851 >             tau_Temp(8) = tau_Temp(8) - d(3) * fy
852 >             tau_Temp(9) = tau_Temp(9) - d(3) * fz
853 >
854               virial_Temp = virial_Temp + &
855                    (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
856  
857            endif
858 <       endif
663 <      
858 >       endif  
859      endif
860  
861 +    
862 +  end subroutine do_eam_pair
863  
667 – end subroutine calc_eam_pair
864  
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 –
865    subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
866  
742 –
867      integer :: atype, nx, j
868      real( kind = DP ), dimension(:) :: xa
869      real( kind = DP ), dimension(:) :: ya
870      real( kind = DP ), dimension(:) :: yppa
871 <    real( kind = DP ) :: x, y, dy, d2y
871 >    real( kind = DP ) :: x, y
872 >    real( kind = DP ) :: dy, d2y
873      real( kind = DP ) :: del, h, a, b, c, d
874 +    integer :: pp_arraySize
875  
876 <
751 <
752 <
753 <    
876 >
877      ! this spline code assumes that the x points are equally spaced
878      ! do not attempt to use this code if they are not.
879      
880      
881      ! find the closest point with a value below our own:
882 <    j = FLOOR(dble(nx-1) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
882 >    j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
883  
884      ! check to make sure we're inside the spline range:
885      if ((j.gt.nx).or.(j.lt.1)) then
886 <       write(default_error,*) "EAM_splint: x is outside bounds of spline"
886 >       write(errMSG,*) "EAM_splint: x is outside bounds of spline"
887 >       call handleError(routineName,errMSG)
888      endif
889      ! check to make sure we haven't screwed up the calculation of j:
890      if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
891         if (j.ne.nx) then
892 <        write(default_error,*) "EAM_splint: x is outside bounding range"
892 >        write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
893 >       call handleError(routineName,errMSG)
894         endif
895      endif
896  
# Line 778 | Line 903 | end subroutine calc_eam_pair
903      d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
904      
905      y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
906 <    
907 <    dy = (ya(j+1)-ya(j))/h &
908 <         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
909 <         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
910 <    
911 <    d2y = a*yppa(j) + b*yppa(j+1)
906 >  
907 >       dy = (ya(j+1)-ya(j))/h &
908 >            - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
909 >            + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
910 >  
911 >  
912 >       d2y = a*yppa(j) + b*yppa(j+1)
913 >  
914  
915    end subroutine eam_splint
916  
917 +
918    subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
919  
792 –  
920  
794 –
921      ! yp1 and ypn are the first derivatives of y at the two endpoints
922      ! if boundary is 'L' the lower derivative is used
923      ! if boundary is 'U' the upper derivative is used
# Line 805 | Line 931 | end subroutine calc_eam_pair
931      real( kind = DP ), dimension(:)        :: yppa
932      real( kind = DP ), dimension(size(xa)) :: u
933      real( kind = DP ) :: yp1,ypn,un,qn,sig,p
934 <    character boundary
934 >    character(len=*) :: boundary
935      
936 <    
936 >    ! make sure the sizes match
937 >    if ((nx /= size(xa)) .or. (nx /= size(ya))) then
938 >       call handleWarning("EAM_SPLINE","Array size mismatch")
939 >    end if
940 >
941      if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
942           (boundary.eq.'b').or.(boundary.eq.'B')) then
943         yppa(1) = -0.5E0_DP

Diff Legend

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