ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/eam.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/eam.F90 (file contents):
Revision 1948 by gezelter, Fri Jan 14 20:31:16 2005 UTC vs.
Revision 2276 by chuckv, Fri Aug 26 20:34:24 2005 UTC

# Line 63 | Line 63 | module eam
63  
64    character(len = 200) :: errMsg
65    character(len=*), parameter :: RoutineName =  "EAM MODULE"
66 < !! Logical that determines if eam arrays should be zeroed
66 >  !! Logical that determines if eam arrays should be zeroed
67    logical :: cleanme = .true.
68    logical :: nmflag  = .false.
69  
70 <
70 >
71    type, private :: EAMtype
72       integer           :: eam_atype      
73       real( kind = DP ) :: eam_dr          
# Line 77 | Line 77 | module eam
77       real( kind = DP ) :: eam_drho      
78       real( kind = DP ) :: eam_rcut    
79       integer           :: eam_atype_map
80 <    
80 >
81       real( kind = DP ), pointer, dimension(:) :: eam_rvals        => null()
82       real( kind = DP ), pointer, dimension(:) :: eam_rhovals      => null()
83       real( kind = DP ), pointer, dimension(:) :: eam_F_rho        => null()
# Line 99 | Line 99 | module eam
99    real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
100  
101  
102 < !! Arrays for MPI storage
102 >  !! Arrays for MPI storage
103   #ifdef IS_MPI
104    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
105    real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row
# Line 115 | Line 115 | module eam
115    type, private :: EAMTypeList
116       integer           :: n_eam_types = 0
117       integer           :: currentAddition = 0
118 <    
118 >
119       type (EAMtype), pointer  :: EAMParams(:) => null()
120 +     integer, pointer         :: atidToEAMType(:) => null()
121    end type EAMTypeList
122  
123  
# Line 132 | Line 133 | module eam
133    public :: calc_eam_prepair_rho
134    public :: calc_eam_preforce_Frho
135    public :: clean_EAM
136 +  public :: destroyEAMTypes
137  
138   contains
139  
# Line 151 | Line 153 | contains
153      integer                                :: eam_ident
154      integer                                :: status
155  
156 <    integer                                :: nAtypes
156 >    integer                                :: nAtypes,nEAMTypes,myATID
157      integer                                :: maxVals
158      integer                                :: alloc_stat
159      integer                                :: current
# Line 162 | Line 164 | contains
164  
165      !! Assume that atypes has already been set and get the total number of types in atypes
166      !! Also assume that every member of atypes is a EAM model.
165 –  
167  
168 +
169      ! check to see if this is the first time into
170      if (.not.associated(EAMList%EAMParams)) then
171 <       call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList)
172 <       EAMList%n_eam_types = nAtypes
173 <       allocate(EAMList%EAMParams(nAtypes))
171 >       call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
172 >       EAMList%n_eam_types = nEAMtypes
173 >       allocate(EAMList%EAMParams(nEAMTypes))
174 >       nAtypes = getSize(atypes)
175 >       allocate(EAMList%atidToEAMType(nAtypes))
176      end if
177  
178      EAMList%currentAddition = EAMList%currentAddition + 1
179      current = EAMList%currentAddition
176 –    
180  
181 +    myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
182 +    EAMList%atidToEAMType(myATID) = current
183 +
184      call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
185      if (alloc_stat /= 0) then
186         status = -1
187         return
188      end if
189  
190 <    ! this is a possible bug, we assume a correspondence between the vector atypes and
185 <    ! EAMAtypes
186 <      
190 >  
191      EAMList%EAMParams(current)%eam_atype    = eam_ident
192      EAMList%EAMParams(current)%eam_lattice  = lattice_constant
193      EAMList%EAMParams(current)%eam_nrho     = eam_nrho
# Line 198 | Line 202 | contains
202    end subroutine newEAMtype
203  
204  
205 +  ! kills all eam types entered and sets simulation to uninitalized
206 +  subroutine destroyEAMtypes()
207 +    integer :: i
208 +    type(EAMType), pointer :: tempEAMType=>null()
209  
210 +    do i = 1, EAMList%n_eam_types
211 +       tempEAMType => eamList%EAMParams(i)
212 +       call deallocate_EAMType(tempEAMType)
213 +    end do
214 +    if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
215 +    eamList%EAMParams => null()
216 +
217 +    eamList%n_eam_types = 0
218 +    eamList%currentAddition = 0
219 +
220 +  end subroutine destroyEAMtypes
221 +
222    subroutine init_EAM_FF(status)
223      integer :: status
224      integer :: i,j
# Line 214 | Line 234 | contains
234         return
235      end if
236  
217 –
218 –       do i = 1, EAMList%currentAddition
237  
238 < ! Build array of r values
238 >    do i = 1, EAMList%currentAddition
239  
240 <          do j = 1,EAMList%EAMParams(i)%eam_nr
223 <             EAMList%EAMParams(i)%eam_rvals(j) = &
224 <                  real(j-1,kind=dp)* &
225 <                  EAMList%EAMParams(i)%eam_dr
226 <              end do
227 < ! Build array of rho values
228 <          do j = 1,EAMList%EAMParams(i)%eam_nrho
229 <             EAMList%EAMParams(i)%eam_rhovals(j) = &
230 <                  real(j-1,kind=dp)* &
231 <                  EAMList%EAMParams(i)%eam_drho
232 <          end do
233 <          ! convert from eV to kcal / mol:
234 <          EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
240 >       ! Build array of r values
241  
242 <          ! precompute the pair potential and get it into kcal / mol:
243 <          EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
244 <          do j = 2, EAMList%EAMParams(i)%eam_nr
245 <             EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
240 <             EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
241 <          enddo
242 >       do j = 1,EAMList%EAMParams(i)%eam_nr
243 >          EAMList%EAMParams(i)%eam_rvals(j) = &
244 >               real(j-1,kind=dp)* &
245 >               EAMList%EAMParams(i)%eam_dr
246         end do
247 <      
247 >       ! Build array of rho values
248 >       do j = 1,EAMList%EAMParams(i)%eam_nrho
249 >          EAMList%EAMParams(i)%eam_rhovals(j) = &
250 >               real(j-1,kind=dp)* &
251 >               EAMList%EAMParams(i)%eam_drho
252 >       end do
253 >       ! convert from eV to kcal / mol:
254 >       EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
255  
256 <       do i = 1,  EAMList%currentAddition
257 <          number_r   = EAMList%EAMParams(i)%eam_nr
258 <          number_rho = EAMList%EAMParams(i)%eam_nrho
259 <          
260 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
250 <               EAMList%EAMParams(i)%eam_rho_r, &
251 <               EAMList%EAMParams(i)%eam_rho_r_pp, &
252 <               0.0E0_DP, 0.0E0_DP, 'N')
253 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
254 <               EAMList%EAMParams(i)%eam_Z_r, &
255 <               EAMList%EAMParams(i)%eam_Z_r_pp, &
256 <               0.0E0_DP, 0.0E0_DP, 'N')
257 <          call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
258 <               EAMList%EAMParams(i)%eam_F_rho, &
259 <               EAMList%EAMParams(i)%eam_F_rho_pp, &
260 <               0.0E0_DP, 0.0E0_DP, 'N')
261 <          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
262 <               EAMList%EAMParams(i)%eam_phi_r, &
263 <               EAMList%EAMParams(i)%eam_phi_r_pp, &
264 <               0.0E0_DP, 0.0E0_DP, 'N')
256 >       ! precompute the pair potential and get it into kcal / mol:
257 >       EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
258 >       do j = 2, EAMList%EAMParams(i)%eam_nr
259 >          EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
260 >          EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
261         enddo
262 +    end do
263  
267 – !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
268 –       !! find the smallest rcut for any eam atype
269 – !       do i = 2, EAMList%currentAddition
270 – !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
271 – !       end do
264  
265 < !       EAM_rcut = current_rcut_max
266 < !       EAM_rcut_orig = current_rcut_max
267 < !       do i = 1, EAMList%currentAddition
268 < !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
269 < !       end do
270 <       !! Allocate arrays for force calculation
271 <      
272 <       call allocateEAM(alloc_stat)
273 <       if (alloc_stat /= 0 ) then
274 <          write(*,*) "allocateEAM failed"
275 <          status = -1
276 <          return
277 <       endif
278 <    
265 >    do i = 1,  EAMList%currentAddition
266 >       number_r   = EAMList%EAMParams(i)%eam_nr
267 >       number_rho = EAMList%EAMParams(i)%eam_nrho
268 >
269 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
270 >            EAMList%EAMParams(i)%eam_rho_r, &
271 >            EAMList%EAMParams(i)%eam_rho_r_pp, &
272 >            0.0E0_DP, 0.0E0_DP, 'N')
273 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
274 >            EAMList%EAMParams(i)%eam_Z_r, &
275 >            EAMList%EAMParams(i)%eam_Z_r_pp, &
276 >            0.0E0_DP, 0.0E0_DP, 'N')
277 >       call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
278 >            EAMList%EAMParams(i)%eam_F_rho, &
279 >            EAMList%EAMParams(i)%eam_F_rho_pp, &
280 >            0.0E0_DP, 0.0E0_DP, 'N')
281 >       call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
282 >            EAMList%EAMParams(i)%eam_phi_r, &
283 >            EAMList%EAMParams(i)%eam_phi_r_pp, &
284 >            0.0E0_DP, 0.0E0_DP, 'N')
285 >    enddo
286 >
287 >    !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
288 >    !! find the smallest rcut for any eam atype
289 >    !       do i = 2, EAMList%currentAddition
290 >    !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
291 >    !       end do
292 >
293 >    !       EAM_rcut = current_rcut_max
294 >    !       EAM_rcut_orig = current_rcut_max
295 >    !       do i = 1, EAMList%currentAddition
296 >    !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
297 >    !       end do
298 >    !! Allocate arrays for force calculation
299 >
300 >    call allocateEAM(alloc_stat)
301 >    if (alloc_stat /= 0 ) then
302 >       write(*,*) "allocateEAM failed"
303 >       status = -1
304 >       return
305 >    endif
306 >
307    end subroutine init_EAM_FF
308  
309 < !! routine checks to see if array is allocated, deallocates array if allocated
310 < !! and then creates the array to the required size
309 >  !! routine checks to see if array is allocated, deallocates array if allocated
310 >  !! and then creates the array to the required size
311    subroutine allocateEAM(status)
312      integer, intent(out) :: status
313  
# Line 330 | Line 350 | contains
350         status = -1
351         return
352      end if
353 <    
353 >
354   #ifdef IS_MPI
355  
356      if (allocated(rho_tmp)) deallocate(rho_tmp)
# Line 367 | Line 387 | contains
387      end if
388  
389  
390 < ! Now do column arrays
390 >    ! Now do column arrays
391  
392      if (allocated(frho_col)) deallocate(frho_col)
393      allocate(frho_col(nAtomsInCol),stat=alloc_stat)
# Line 393 | Line 413 | contains
413         status = -1
414         return
415      end if
416 <  
416 >
417   #endif
418  
419    end subroutine allocateEAM
420  
421 < !! C sets rcut to be the largest cutoff of any atype
422 < !! present in this simulation. Doesn't include all atypes
423 < !! sim knows about, just those in the simulation.
421 >  !! C sets rcut to be the largest cutoff of any atype
422 >  !! present in this simulation. Doesn't include all atypes
423 >  !! sim knows about, just those in the simulation.
424    subroutine setCutoffEAM(rcut, status)
425      real(kind=dp) :: rcut
426      integer :: status
# Line 413 | Line 433 | contains
433  
434  
435    subroutine clean_EAM()
436 <  
437 <   ! clean non-IS_MPI first
436 >
437 >    ! clean non-IS_MPI first
438      frho = 0.0_dp
439      rho  = 0.0_dp
440      dfrhodrho = 0.0_dp
441 < ! clean MPI if needed
441 >    ! clean MPI if needed
442   #ifdef IS_MPI
443      frho_row = 0.0_dp
444      frho_col = 0.0_dp
# Line 442 | Line 462 | contains
462  
463  
464      if (present(stat)) stat = 0
465 <    
465 >
466      allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
467      if (alloc_stat /= 0 ) then
468         if (present(stat)) stat = -1
# Line 493 | Line 513 | contains
513         if (present(stat)) stat = -1
514         return
515      end if
496 –      
516  
517 +
518    end subroutine allocate_EAMType
519  
520  
# Line 502 | Line 522 | contains
522      type (EAMtype), pointer :: thisEAMType
523  
524      ! free Arrays in reverse order of allocation...
525 <    deallocate(thisEAMType%eam_phi_r_pp)      
526 <    deallocate(thisEAMType%eam_rho_r_pp)  
527 <    deallocate(thisEAMType%eam_Z_r_pp)  
528 <    deallocate(thisEAMType%eam_F_rho_pp)  
529 <    deallocate(thisEAMType%eam_phi_r)      
530 <    deallocate(thisEAMType%eam_rho_r)      
531 <    deallocate(thisEAMType%eam_Z_r)  
532 <    deallocate(thisEAMType%eam_F_rho)
533 <    deallocate(thisEAMType%eam_rhovals)
534 <    deallocate(thisEAMType%eam_rvals)
535 <  
525 >    if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)      
526 >    if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)  
527 >    if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)  
528 >    if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)  
529 >    if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)      
530 >    if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)      
531 >    if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)  
532 >    if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
533 >    if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
534 >    if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
535 >
536    end subroutine deallocate_EAMType
537  
538 < !! Calculates rho_r
538 >  !! Calculates rho_r
539    subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
540      integer :: atom1,atom2
541      real(kind = dp), dimension(3) :: d
# Line 529 | Line 549 | contains
549      ! we don't use the derivatives, dummy variables
550      real( kind = dp) :: drho,d2rho
551      integer :: eam_err
552 <  
552 >
553      integer :: myid_atom1
554      integer :: myid_atom2
555  
556 < ! check to see if we need to be cleaned at the start of a force loop
537 <      
556 >    ! check to see if we need to be cleaned at the start of a force loop
557  
558  
559  
560 +
561   #ifdef IS_MPI
562      myid_atom1 = atid_Row(atom1)
563      myid_atom2 = atid_Col(atom2)
# Line 557 | Line 577 | contains
577              r, rho_i_at_j,drho,d2rho)
578  
579  
580 <      
580 >
581   #ifdef  IS_MPI
582         rho_col(atom2) = rho_col(atom2) + rho_i_at_j
583   #else
584         rho(atom2) = rho(atom2) + rho_i_at_j
585   #endif
586 < !       write(*,*) atom1,atom2,r,rho_i_at_j
587 <       endif
586 >       !       write(*,*) atom1,atom2,r,rho_i_at_j
587 >    endif
588  
589 <       if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
590 <          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
591 <               EAMList%EAMParams(myid_atom2)%eam_rvals, &
592 <               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
593 <               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
594 <               r, rho_j_at_i,drho,d2rho)
589 >    if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
590 >       call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
591 >            EAMList%EAMParams(myid_atom2)%eam_rvals, &
592 >            EAMList%EAMParams(myid_atom2)%eam_rho_r, &
593 >            EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
594 >            r, rho_j_at_i,drho,d2rho)
595  
596 <    
597 <      
598 <      
596 >
597 >
598 >
599   #ifdef  IS_MPI
600 <          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
600 >       rho_row(atom1) = rho_row(atom1) + rho_j_at_i
601   #else
602 <          rho(atom1) = rho(atom1) + rho_j_at_i
602 >       rho(atom1) = rho(atom1) + rho_j_at_i
603   #endif
604 <       endif
604 >    endif
605  
606  
607  
# Line 604 | Line 624 | contains
624      integer :: me
625      integer :: n_rho_points
626  
627 <  
627 >
628      cleanme = .true.
629      !! Scatter the electron density from  pre-pair calculation back to local atoms
630   #ifdef IS_MPI
631      call scatter(rho_row,rho,plan_atom_row,eam_err)
632      if (eam_err /= 0 ) then
633 <      write(errMsg,*) " Error scattering rho_row into rho"
634 <      call handleError(RoutineName,errMesg)
635 <   endif      
633 >       write(errMsg,*) " Error scattering rho_row into rho"
634 >       call handleError(RoutineName,errMesg)
635 >    endif
636      call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
637      if (eam_err /= 0 ) then
638 <      write(errMsg,*) " Error scattering rho_col into rho"
639 <      call handleError(RoutineName,errMesg)
640 <   endif
638 >       write(errMsg,*) " Error scattering rho_col into rho"
639 >       call handleError(RoutineName,errMesg)
640 >    endif
641  
642 <      rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
642 >    rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
643   #endif
644  
645  
646  
647 < !! Calculate F(rho) and derivative
647 >    !! Calculate F(rho) and derivative
648      do atom = 1, nlocal
649         me = atid(atom)
650         n_rho_points = EAMList%EAMParams(me)%eam_nrho
# Line 654 | Line 674 | contains
674  
675      enddo
676  
657 –  
677  
678 +
679   #ifdef IS_MPI
680      !! communicate f(rho) and derivatives back into row and column arrays
681      call gather(frho,frho_row,plan_atom_row, eam_err)
# Line 685 | Line 705 | contains
705      endif
706   #endif
707  
708 <  
708 >
709    end subroutine calc_eam_preforce_Frho
710  
711  
# Line 703 | Line 723 | contains
723      real( kind = dp ), intent(inout), dimension(3) :: fpair
724  
725      logical, intent(in) :: do_pot
726 <    
726 >
727      real( kind = dp ) :: drdx,drdy,drdz
728      real( kind = dp ) :: d2
729      real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
# Line 721 | Line 741 | contains
741      integer  :: mytype_atom1
742      integer  :: mytype_atom2
743  
744 < !Local Variables
725 <    
726 <   ! write(*,*) "Frho: ", Frho(atom1)
744 >    !Local Variables
745  
746 +    ! write(*,*) "Frho: ", Frho(atom1)
747 +
748      phab = 0.0E0_DP
749      dvpdr = 0.0E0_DP
750      d2vpdrdr = 0.0E0_DP
# Line 742 | Line 762 | contains
762         rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
763         ! get type specific cutoff for atom 2
764         rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
765 <      
765 >
766         drdx = d(1)/rij
767         drdy = d(2)/rij
768         drdz = d(3)/rij
769 <      
769 >
770         if (rij.lt.rci) then
771            call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
772                 EAMList%EAMParams(mytype_atom1)%eam_rvals, &
# Line 768 | Line 788 | contains
788                 EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
789                 EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
790                 rij, rhb,drhb,d2rhb)
791 <          
791 >
792            !! Calculate Phi(r) for atom2.
793            call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
794                 EAMList%EAMParams(mytype_atom2)%eam_rvals, &
# Line 786 | Line 806 | contains
806                 pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
807                 (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
808         endif
809 <      
809 >
810         if (rij.lt.rcj) then
811            phab = phab + 0.5E0_DP*(rha/rhb)*phb
812            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
# Line 796 | Line 816 | contains
816                 phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
817                 (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
818         endif
819 <      
819 >
820         drhoidr = drha
821         drhojdr = drhb
822  
# Line 811 | Line 831 | contains
831   #else
832         dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
833              + dvpdr
834 <      ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
834 >       ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
835   #endif
836  
837         fx = dudr * drdx
# Line 828 | Line 848 | contains
848         f_Row(1,atom1) = f_Row(1,atom1) + fx
849         f_Row(2,atom1) = f_Row(2,atom1) + fy
850         f_Row(3,atom1) = f_Row(3,atom1) + fz
851 <      
851 >
852         f_Col(1,atom2) = f_Col(1,atom2) - fx
853         f_Col(2,atom2) = f_Col(2,atom2) - fy
854         f_Col(3,atom2) = f_Col(3,atom2) - fz
# Line 841 | Line 861 | contains
861         f(1,atom1) = f(1,atom1) + fx
862         f(2,atom1) = f(2,atom1) + fy
863         f(3,atom1) = f(3,atom1) + fz
864 <      
864 >
865         f(1,atom2) = f(1,atom2) - fx
866         f(2,atom2) = f(2,atom2) - fy
867         f(3,atom2) = f(3,atom2) - fz
868   #endif
869 <      
869 >
870         vpair = vpair + phab
871   #ifdef IS_MPI
872         id1 = AtomRowToGlobal(atom1)
# Line 855 | Line 875 | contains
875         id1 = atom1
876         id2 = atom2
877   #endif
878 <      
878 >
879         if (molMembershipList(id1) .ne. molMembershipList(id2)) then
880 <          
880 >
881            fpair(1) = fpair(1) + fx
882            fpair(2) = fpair(2) + fy
883            fpair(3) = fpair(3) + fz
884 <          
884 >
885         endif
886 <    
886 >
887         if (nmflag) then
888  
889            drhoidr = drha
# Line 877 | Line 897 | contains
897                 d2rhojdrdr*dfrhodrho_row(atom1) + &
898                 drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
899                 drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
900 <              
900 >
901   #else
902  
903            d2 = d2vpdrdr + &
# Line 887 | Line 907 | contains
907                 drhojdr*drhojdr*d2frhodrhodrho(atom1)
908   #endif
909         end if
910 <      
911 <    endif    
910 >
911 >    endif
912    end subroutine do_eam_pair
913  
914  
# Line 903 | Line 923 | contains
923      real( kind = DP ) :: del, h, a, b, c, d
924      integer :: pp_arraySize
925  
926 <
926 >
927      ! this spline code assumes that the x points are equally spaced
928      ! do not attempt to use this code if they are not.
929 <    
930 <    
929 >
930 >
931      ! find the closest point with a value below our own:
932      j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
933  
# Line 919 | Line 939 | contains
939      ! check to make sure we haven't screwed up the calculation of j:
940      if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
941         if (j.ne.nx) then
942 <        write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
943 <       call handleError(routineName,errMSG)
942 >          write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
943 >          call handleError(routineName,errMSG)
944         endif
945      endif
946  
947      del = xa(j+1) - x
948      h = xa(j+1) - xa(j)
949 <    
949 >
950      a = del / h
951      b = 1.0E0_DP - a
952      c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
953      d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
954 <    
954 >
955      y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
936 –  
937 –       dy = (ya(j+1)-ya(j))/h &
938 –            - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
939 –            + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
940 –  
941 –  
942 –       d2y = a*yppa(j) + b*yppa(j+1)
943 –  
956  
957 +    dy = (ya(j+1)-ya(j))/h &
958 +         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
959 +         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
960 +
961 +
962 +    d2y = a*yppa(j) + b*yppa(j+1)
963 +
964 +
965    end subroutine eam_splint
966  
967  
# Line 953 | Line 973 | contains
973      ! if boundary is 'U' the upper derivative is used
974      ! if boundary is 'B' then both derivatives are used
975      ! if boundary is anything else, then both derivatives are assumed to be 0
976 <    
976 >
977      integer :: nx, i, k, max_array_size
978 <    
978 >
979      real( kind = DP ), dimension(:)        :: xa
980      real( kind = DP ), dimension(:)        :: ya
981      real( kind = DP ), dimension(:)        :: yppa
982      real( kind = DP ), dimension(size(xa)) :: u
983      real( kind = DP ) :: yp1,ypn,un,qn,sig,p
984      character(len=*) :: boundary
985 <    
985 >
986      ! make sure the sizes match
987      if ((nx /= size(xa)) .or. (nx /= size(ya))) then
988         call handleWarning("EAM_SPLINE","Array size mismatch")
# Line 977 | Line 997 | contains
997         yppa(1) = 0.0E0_DP
998         u(1)  = 0.0E0_DP
999      endif
1000 <    
1000 >
1001      do i = 2, nx - 1
1002         sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1003         p = sig * yppa(i-1) + 2.0E0_DP
# Line 986 | Line 1006 | contains
1006              (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1007              (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1008      enddo
1009 <    
1009 >
1010      if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1011           (boundary.eq.'b').or.(boundary.eq.'B')) then
1012         qn = 0.5E0_DP
# Line 998 | Line 1018 | contains
1018      endif
1019  
1020      yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1021 <    
1021 >
1022      do k = nx-1, 1, -1
1023         yppa(k)=yppa(k)*yppa(k+1)+u(k)
1024      enddo

Diff Legend

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