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 492 by chuckv, Mon Apr 14 18:19:15 2003 UTC vs.
Revision 631 by chuckv, Thu Jul 17 19:25:51 2003 UTC

# Line 1 | Line 1
1 < module calc_eam
2 <  use definitions, ONLY : DP
1 > module eam
2 >  use definitions, ONLY : DP,default_error
3 >  use simulation
4    use force_globals
5 +  use status
6 +  use atype_module
7   #ifdef MPI
8    use mpiSimulation
9   #endif
10 +  implicit none
11 +  PRIVATE
12  
13  
14 +  logical, save :: EAM_FF_initialized = .false.
15 +  integer, save :: EAM_Mixing_Policy
16 +  real(kind = dp), save :: EAM_rcut
17  
18 +  character(len = 200) :: errMsg
19 +  character(len=*), parameter :: RoutineName =  "EAM MODULE"
20 +  logical :: cleanme = .true.
21  
11 –  !! standard eam stuff  
12 –  integer                                        :: n_eam_atypes
13 –  integer,           allocatable, dimension(:)   :: eam_atype    
14 –  real( kind = DP ), allocatable, dimension(:)   :: eam_dr      
15 –  integer, allocatable, dimension(:)   :: eam_nr      
16 –  integer, allocatable, dimension(:)   :: eam_nrho      
17 –  real( kind = DP ), allocatable, dimension(:)   :: eam_lattice      
18 –  real( kind = DP ), allocatable, dimension(:)   :: eam_drho    
19 –  integer          , allocatable, dimension(:)   :: eam_atype_map
20 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_rvals    
21 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_rhovals  
22 –  real( kind = DP ), allocatable, dimension(:)   :: eam_rcut    
23 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_F_rho    
24 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_Z_r      
25 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_rho_r    
26 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_phi_r    
27 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_F_rho_pp
28 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_Z_r_pp  
29 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_rho_r_pp
30 –  real( kind = DP ), allocatable, dimension(:,:) :: eam_phi_r_pp
22  
23 +
24 +  type, private :: EAMtype
25 +     integer           :: eam_atype      
26 +     real( kind = DP ) :: eam_dr          
27 +     integer           :: eam_nr          
28 +     integer           :: eam_nrho          
29 +     real( kind = DP ) :: eam_lattice        
30 +     real( kind = DP ) :: eam_drho      
31 +     real( kind = DP ) :: eam_rcut    
32 +     integer           :: eam_atype_map
33 +    
34 +     real( kind = DP ), pointer, dimension(:) :: eam_rvals        => null()
35 +     real( kind = DP ), pointer, dimension(:) :: eam_rhovals      => null()
36 +     real( kind = DP ), pointer, dimension(:) :: eam_F_rho        => null()
37 +     real( kind = DP ), pointer, dimension(:) :: eam_Z_r          => null()
38 +     real( kind = DP ), pointer, dimension(:) :: eam_rho_r        => null()
39 +     real( kind = DP ), pointer, dimension(:) :: eam_phi_r        => null()
40 +     real( kind = DP ), pointer, dimension(:) :: eam_F_rho_pp     => null()
41 +     real( kind = DP ), pointer, dimension(:) :: eam_Z_r_pp       => null()
42 +     real( kind = DP ), pointer, dimension(:) :: eam_rho_r_pp     => null()
43 +     real( kind = DP ), pointer, dimension(:) :: eam_phi_r_pp     => null()
44 +  end type EAMtype
45  
33 –  real( kind = DP ), private :: time0,time1,time2,time3
46  
47 <  integer, private :: eam_err
47 >  !! Arrays for derivatives used in force calculation
48 >  real( kind = dp), dimension(:), allocatable :: frho
49 >  real( kind = dp), dimension(:), allocatable :: rho
50  
51 <  private :: mass_weight
52 <  private :: allocate_eam_atype,allocate_eam_module,deallocate_eam_module
39 <  private :: read_eam_pot, get_eam_sizes
51 >  real( kind = dp), dimension(:), allocatable :: dfrhodrho
52 > !  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
53  
54  
55 < contains
56 <  subroutine initialize_eam()
55 > !! Arrays for MPI storage
56 > #ifdef MPI
57 >  real( kind = dp), dimension(:), allocatable :: dfrhodrho_col
58 >  real( kind = dp), dimension(:), allocatable :: dfrhodrho_row
59 >  real( kind = dp), dimension(:), allocatable :: frho_row
60 >  real( kind = dp), dimension(:), allocatable :: frho_col
61 >  real( kind = dp), dimension(:), allocatable :: rho_row
62 >  real( kind = dp), dimension(:), allocatable :: rho_col
63  
64 + #endif
65  
66 +  type, private :: EAMTypeList
67 +     integer           :: n_eam_types = 0
68 +     integer           :: currentAddition = 0
69 +    
70 +     type (EAMtype), pointer  :: EAMParams(:) => null()
71 +  end type EAMTypeList
72  
47 –    character(len=80) :: eam_pot_file
48 –    integer :: i, j,  max_size, prev_max_size
49 –    integer :: number_rho, number_r
50 –    integer :: eam_unit
51 –    integer :: this_error
52 –    character(len=300) :: msg
53 –    integer, external :: nfiles
54 –    !for mpi
73  
74 +  type (eamTypeList) :: EAMList
75  
76 < #ifdef MPI
58 <    if (node == 0)  &  
59 <         n_eam_atypes = nfiles(trim(eam_pot_dir)//char(0))  
76 >  !! standard eam stuff  
77  
61 –    call mpi_bcast(n_eam_atypes,1,mpi_integer,0,mpi_comm_world,mpi_err)
62 –    if (n_eam_atypes == -1) then
63 –       call error("INITIALIZE_EAM","NO EAM potentials found!")
64 –    endif
65 –    write(msg,'(a5,i4,a12,i5,a14)') 'Node: ',node,' reading ...', &
66 –         n_eam_atypes, ' eam atom types'
67 –    call info('INITIALIZE_EAM', trim(msg))
68 – #else
69 –    n_eam_atypes = nfiles(trim(eam_pot_dir)//char(0))
70 –    if (n_eam_atypes == -1) then
71 –       call error("INITIALIZE_EAM","NO EAM potentials found!")
72 –    endif
78  
79 <    write(msg,'(a12,i5,a14)') ' Reading ...', &
80 <         n_eam_atypes, ' eam atom types'
81 <    call info('INITIALIZE_EAM', trim(msg))
82 < #endif
79 >  public :: init_EAM_FF
80 > !  public :: EAM_new_rcut
81 > !  public :: do_EAM_pair
82 >  public :: newEAMtype
83  
84 +  
85  
86 <    call allocate_eam_atype(n_eam_atypes)
86 > contains
87  
88  
89 +  subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
90 +       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
91 +       eam_ident,status)
92 +    real (kind = dp )                      :: lattice_constant
93 +    integer                                :: eam_nrho
94 +    real (kind = dp )                      :: eam_drho
95 +    integer                                :: eam_nr
96 +    real (kind = dp )                      :: eam_dr
97 +    real (kind = dp )                      :: rcut
98 +    real (kind = dp ), dimension(eam_nr)   :: eam_Z_r
99 +    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
100 +    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
101 +    integer                                :: eam_ident
102 +    integer                                :: status
103  
104 <    !! get largest number of data points for any potential
105 < #ifdef MPI
106 <    if (node == 0) then
107 < #endif
108 <       prev_max_size = 0
109 <       do i = 1, n_eam_atypes
110 <          call getfilename(i, eam_pot_file)
111 <          max_size = max(get_eam_sizes( &
112 <               trim(eam_pot_dir) // '/' // eam_pot_file), &
113 <               prev_max_size)
114 <          prev_max_size = max_size
115 <       end do
116 < #ifdef MPI
104 >    integer                                :: nAtypes
105 >    integer                                :: maxVals
106 >    integer                                :: alloc_stat
107 >    integer                                :: current
108 >    integer,pointer                        :: Matchlist(:) => null()
109 >    status = 0
110 >
111 >    !! Assume that atypes has already been set and get the total number of types in atypes
112 >    !! Also assume that every member of atypes is a EAM model.
113 >  
114 >
115 >    ! check to see if this is the first time into
116 >    if (.not.associated(EAMList%EAMParams)) then
117 >       call getMatchingElementList(atypes, "is_EAM", .true., nAtypes, MatchList)
118 >       EAMList%n_eam_types = nAtypes
119 >       allocate(EAMList%EAMParams(nAtypes))
120      end if
121  
122 +    EAMList%currentAddition = EAMList%currentAddition + 1
123 +    current = EAMList%currentAddition
124 +    
125  
126 <    call mpi_bcast(max_size,1,mpi_integer,0,mpi_comm_world,mpi_err)
127 < #endif
126 >    !call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
127 >    if (alloc_stat /= 0) then
128 >       status = -1
129 >       return
130 >    end if
131  
103 –    call allocate_eam_module(n_eam_atypes,max_size)
104 –    allocate(eam_atype_map(get_max_atype()))
132  
133 < #ifdef MPI
134 <    if (node == 0) then
135 < #endif
136 <       do i = 1, n_eam_atypes
137 <          call getfilename(i, eam_pot_file)
138 <          call read_eam_pot(i,trim(eam_pot_dir) // '/' // eam_pot_file, &
139 <               this_error)
133 >    ! this is a possible bug, we assume a correspondence between the vector atypes and
134 >    ! EAMAtypes
135 >      
136 >    EAMList%EAMParams(current)%eam_atype    = eam_ident
137 >    EAMList%EAMParams(current)%eam_lattice  = lattice_constant
138 >    EAMList%EAMParams(current)%eam_nrho     = eam_nrho
139 >    EAMList%EAMParams(current)%eam_drho     = eam_drho
140 >    EAMList%EAMParams(current)%eam_nr       = eam_nr
141 >    EAMList%EAMParams(current)%eam_rcut     = rcut
142 >    EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
143 >    EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
144 >    EAMList%EAMParams(current)%eam_F_rho    = eam_F_rho
145  
146 <          do j = 1, eam_nr(i)
115 <             eam_rvals(j,i) = dble(j-1)*eam_dr(i)
116 <          enddo
146 >  end subroutine newEAMtype
147  
118 –          do j = 1, eam_nrho(i)
119 –             eam_rhovals(j,i) = dble(j-1)*eam_drho(i)
120 –          enddo
148  
149 +
150 +  subroutine init_EAM_FF(status)
151 +    integer :: status
152 +    integer :: i,j
153 +    real(kind=dp) :: current_rcut_max
154 +    integer :: alloc_stat
155 +    integer :: number_r, number_rho
156 +
157 +       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
162 +
163 +          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
166 +
167            ! convert from eV to kcal / mol:
168 <          do j = 1, eam_nrho(i)
124 <             eam_F_rho(j,i) = eam_F_rho(j,i)*23.06054E0_DP
125 <          enddo
168 >          EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
169  
170            ! precompute the pair potential and get it into kcal / mol:
171 <          eam_phi_r(1,i) = 0.0E0_DP
172 <          do j = 2, eam_nr(i)
173 <             eam_phi_r(j,i) = (eam_Z_r(j,i)**2)/eam_rvals(j,i)
174 <             eam_phi_r(j,i) = eam_phi_r(j,i)*331.999296E0_DP
171 >          EAMList%EAMParams(i)%eam_phi_r = 0.0E0_DP
172 >          do j = 2, EAMList%EAMParams(i)%eam_nr
173 >             EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
174 >             EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
175            enddo
176 +       end do
177  
178 +       do i = 1,  EAMList%currentAddition
179 +          number_r   = EAMList%EAMParams(i)%eam_nr
180 +          number_rho = EAMList%EAMParams(i)%eam_nrho
181 +          
182 +          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
183 +               EAMList%EAMParams(i)%eam_rho_r, &
184 +               EAMList%EAMParams(i)%eam_rho_r_pp, &
185 +               0.0E0_DP, 0.0E0_DP, 'N')
186 +          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
187 +               EAMList%EAMParams(i)%eam_Z_r, &
188 +               EAMList%EAMParams(i)%eam_Z_r_pp, &
189 +               0.0E0_DP, 0.0E0_DP, 'N')
190 +          call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
191 +               EAMList%EAMParams(i)%eam_F_rho, &
192 +               EAMList%EAMParams(i)%eam_F_rho_pp, &
193 +               0.0E0_DP, 0.0E0_DP, 'N')
194 +          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
195 +               EAMList%EAMParams(i)%eam_phi_r, &
196 +               EAMList%EAMParams(i)%eam_phi_r_pp, &
197 +               0.0E0_DP, 0.0E0_DP, 'N')
198 +       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
135 – #ifdef MPI
136 –       call info('INITIALIZE_EAM','NODE 0: Distributing spline arrays')
137 –    endif
205  
206 <    call mpi_bcast(this_error,n_eam_atypes,mpi_integer,0, &
207 <         mpi_comm_world,mpi_err)
208 <    if (this_error /= 0) then
209 <       call error('INITIALIZE_EAM',"Cannot read eam files")
143 <    endif
206 >       EAM_rcut = current_rcut_max
207 > !       do i = 1, EAMList%currentAddition
208 > !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
209 > !       end do
210  
145 –    call mpi_bcast(eam_atype,n_eam_atypes,mpi_integer,0, &
146 –         mpi_comm_world,mpi_err)
211  
212 <    !! distribute values to cluster......
213 <    call mpi_bcast(eam_nr,n_eam_atypes,mpi_integer,&
214 <         0,mpi_comm_world,mpi_err)
215 <    call mpi_bcast(eam_nrho,n_eam_atypes,mpi_integer,&
216 <         0,mpi_comm_world,mpi_err)
217 <    call mpi_bcast(eam_rvals,n_eam_atypes*max_size,mpi_double_precision, &
154 <         0,mpi_comm_world,mpi_err)
155 <    call mpi_bcast(eam_rcut,n_eam_atypes,mpi_double_precision, &
156 <         0,mpi_comm_world,mpi_err)
157 <    call mpi_bcast(eam_rhovals,n_eam_atypes*max_size,mpi_double_precision, &
158 <         0,mpi_comm_world,mpi_err)
212 >       !! Allocate arrays for force calculation
213 >          call allocateEAM(alloc_stat)
214 >          if (alloc_stat /= 0 ) then
215 >             status = -1
216 >             return
217 >          endif
218  
219 <    !! distribute arrays
161 <    call mpi_bcast(eam_rho_r,n_eam_atypes*max_size,mpi_double_precision, &
162 <         0,mpi_comm_world,mpi_err)
163 <    call mpi_bcast(eam_Z_r,n_eam_atypes*max_size,mpi_double_precision, &
164 <         0,mpi_comm_world,mpi_err)
165 <    call mpi_bcast(eam_F_rho,n_eam_atypes*max_size,mpi_double_precision, &
166 <         0,mpi_comm_world,mpi_err)
167 <    call mpi_bcast(eam_phi_r,n_eam_atypes*max_size,mpi_double_precision, &
168 <         0,mpi_comm_world,mpi_err)
219 >  end subroutine init_EAM_FF
220  
221 + !! routine checks to see if array is allocated, deallocates array if allocated
222 + !! and then creates the array to the required size
223 +  subroutine allocateEAM(status)
224 +    integer, intent(out) :: status
225 +
226 +    integer :: nlocal
227 + #ifdef MPI
228 +    integer :: nrow
229 +    integer :: ncol
230   #endif
231 <    call info('INITIALIZE_EAM', 'creating splines')
231 >    integer :: alloc_stat
232  
173 –    do i = 1, n_eam_atypes
174 –       number_r = eam_nr(i)
175 –       number_rho = eam_nrho(i)
233  
234 <       call eam_spline(i, number_r, eam_rvals, eam_rho_r, eam_rho_r_pp, &
178 <            0.0E0_DP, 0.0E0_DP, 'N')
179 <       call eam_spline(i, number_r, eam_rvals, eam_Z_r, eam_Z_r_pp, &
180 <            0.0E0_DP, 0.0E0_DP, 'N')
181 <       call eam_spline(i, number_rho, eam_rhovals, eam_F_rho, eam_F_rho_pp, &
182 <            0.0E0_DP, 0.0E0_DP, 'N')
183 <       call eam_spline(i, number_r, eam_rvals, eam_phi_r, eam_phi_r_pp, &
184 <            0.0E0_DP, 0.0E0_DP, 'N')
185 <    enddo
234 >    nlocal = getNlocal()
235  
236 <    do i = 1, n_eam_atypes
237 <       eam_atype_map(eam_atype(i)) = i
238 <    end do
236 > #ifdef MPI
237 >    nrow = getNrow(plan_row)
238 >    ncol = getNcol(plan_col)
239 > #endif
240  
241 +    if (allocated(frho)) deallocate(frho)
242 +    allocate(frho(nlocal),stat=alloc_stat)
243 +    if (alloc_stat /= 0) then
244 +       status = -1
245 +       return
246 +    end if
247 +    if (allocated(rho)) deallocate(rho)
248 +    allocate(rho(nlocal),stat=alloc_stat)
249 +    if (alloc_stat /= 0) then
250 +       status = -1
251 +       return
252 +    end if
253 +    if (allocated(dfrhodrho)) deallocate(dfrhodrho)
254 +    allocate(dfrhodrho(nlocal),stat=alloc_stat)
255 +    if (alloc_stat /= 0) then
256 +       status = -1
257 +       return
258 +    end if
259 +    
260 + #ifdef MPI
261  
262 +    if (allocated(frho_row)) deallocate(frho_row)
263 +    allocate(frho_row(nrow),stat=alloc_stat)
264 +    if (alloc_stat /= 0) then
265 +       status = -1
266 +       return
267 +    end if
268 +    if (allocated(rho_row)) deallocate(rho_row)
269 +    allocate(rho_row(nrow),stat=alloc_stat)
270 +    if (alloc_stat /= 0) then
271 +       status = -1
272 +       return
273 +    end if
274 +    if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
275 +    allocate(dfrhodrho_row(nrow),stat=alloc_stat)
276 +    if (alloc_stat /= 0) then
277 +       status = -1
278 +       return
279 +    end if
280  
281 <    call info('INITIALIZE_EAM','Done creating splines')
281 > ! Now do column arrays
282  
283 <    return
284 <  end subroutine initialize_eam
283 >    if (allocated(frho_col)) deallocate(frho_col)
284 >    allocate(frho_col(ncol),stat=alloc_stat)
285 >    if (alloc_stat /= 0) then
286 >       status = -1
287 >       return
288 >    end if
289 >    if (allocated(rho_col)) deallocate(rho_col)
290 >    allocate(rho_col(ncol),stat=alloc_stat)
291 >    if (alloc_stat /= 0) then
292 >       status = -1
293 >       return
294 >    end if
295 >    if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
296 >    allocate(dfrhodrho_col(ncol),stat=alloc_stat)
297 >    if (alloc_stat /= 0) then
298 >       status = -1
299 >       return
300 >    end if
301  
302 + #endif
303  
304 <  subroutine allocate_eam_atype(n_size_atype)
200 <    integer, intent(in) :: n_size_atype
304 >  end subroutine allocateEAM
305  
202 –    allocate(eam_atype(n_size_atype))  
203 –    allocate(eam_drho(n_size_atype))    
204 –    allocate(eam_dr(n_size_atype))      
205 –    allocate(eam_nr(n_size_atype))
206 –    allocate(eam_nrho(n_size_atype))
207 –    allocate(eam_lattice(n_size_atype))      
208 –    allocate(eam_rcut(n_size_atype))  
306  
307 <  end subroutine allocate_eam_atype
307 >  subroutine clean_EAM()
308  
309 <  subroutine allocate_eam_module(n_size_atype,n_eam_points)
310 <    integer, intent(in) :: n_eam_points
311 <    integer, intent(in) :: n_size_atype
309 > ! clean non-MPI first
310 >    frho = 0.0_dp
311 >    rho  = 0.0_dp
312 >    dfrhodrho = 0.0_dp
313 > ! clean MPI if needed
314 > #ifdef MPI
315 >    frho_row = 0.0_dp
316 >    frho_col = 0.0_dp
317 >    rho_row  = 0.0_dp
318 >    rho_col  = 0.0_dp
319 >    dfrhodrho_row = 0.0_dp
320 >    dfrhodrho_col = 0.0_dp
321 > #endif
322 >  end subroutine clean_EAM
323  
216 –    allocate(eam_rvals(n_eam_points,n_size_atype))  
217 –    allocate(eam_rhovals(n_eam_points,n_size_atype))  
218 –    allocate(eam_F_rho(n_eam_points,n_size_atype))  
219 –    allocate(eam_Z_r(n_eam_points,n_size_atype))        
220 –    allocate(eam_rho_r(n_eam_points,n_size_atype))      
221 –    allocate(eam_phi_r(n_eam_points,n_size_atype))      
222 –    allocate(eam_F_rho_pp(n_eam_points,n_size_atype))  
223 –    allocate(eam_Z_r_pp(n_eam_points,n_size_atype))  
224 –    allocate(eam_rho_r_pp(n_eam_points,n_size_atype))  
225 –    allocate(eam_phi_r_pp(n_eam_points,n_size_atype))  
324  
227 –  end subroutine allocate_eam_module
325  
326 <  subroutine deallocate_eam_module()
326 >  subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
327 >    integer, intent(in)          :: eam_n_rho
328 >    integer, intent(in)          :: eam_n_r
329 >    type (EAMType), pointer      :: thisEAMType
330 >    integer, optional   :: stat
331 >    integer             :: alloc_stat
332  
231 –    deallocate(eam_atype)  
232 –    deallocate(eam_drho)    
233 –    deallocate(eam_dr)      
234 –    deallocate(eam_nr)
235 –    deallocate(eam_nrho)
236 –    deallocate(eam_lattice)      
237 –    deallocate(eam_atype_map)
238 –    deallocate(eam_rvals)  
239 –    deallocate(eam_rhovals)  
240 –    deallocate(eam_rcut)  
241 –    deallocate(eam_Z_r)        
242 –    deallocate(eam_rho_r)      
243 –    deallocate(eam_phi_r)      
244 –    deallocate(eam_F_rho_pp)  
245 –    deallocate(eam_Z_r_pp)  
246 –    deallocate(eam_rho_r_pp)  
247 –    deallocate(eam_phi_r_pp)  
333  
334 <  end subroutine deallocate_eam_module
334 >
335 >    if (present(stat)) stat = 0
336 >    
337 >    allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)  
338 >    if (alloc_stat /= 0 ) then
339 >       if (present(stat)) stat = -1
340 >       return
341 >    end if
342 >    allocate(thisEAMType%eam_rhovals(eam_n_rho),stat=alloc_stat)  
343 >    if (alloc_stat /= 0 ) then
344 >       if (present(stat)) stat = -1
345 >       return
346 >    end if
347 >    allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat)  
348 >    if (alloc_stat /= 0 ) then
349 >       if (present(stat)) stat = -1
350 >       return
351 >    end if
352 >    allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat)        
353 >    if (alloc_stat /= 0 ) then
354 >       if (present(stat)) stat = -1
355 >       return
356 >    end if
357 >    allocate(thisEAMType%eam_rho_r(eam_n_r),stat=alloc_stat)      
358 >    if (alloc_stat /= 0 ) then
359 >       if (present(stat)) stat = -1
360 >       return
361 >    end if
362 >    allocate(thisEAMType%eam_phi_r(eam_n_r),stat=alloc_stat)      
363 >    if (alloc_stat /= 0 ) then
364 >       if (present(stat)) stat = -1
365 >       return
366 >    end if
367 >    allocate(thisEAMType%eam_F_rho_pp(eam_n_rho),stat=alloc_stat)  
368 >    if (alloc_stat /= 0 ) then
369 >       if (present(stat)) stat = -1
370 >       return
371 >    end if
372 >    allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat)  
373 >    if (alloc_stat /= 0 ) then
374 >       if (present(stat)) stat = -1
375 >       return
376 >    end if
377 >    allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat)  
378 >    if (alloc_stat /= 0 ) then
379 >       if (present(stat)) stat = -1
380 >       return
381 >    end if
382 >    allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat)
383 >    if (alloc_stat /= 0 ) then
384 >       if (present(stat)) stat = -1
385 >       return
386 >    end if
387 >      
388 >
389 >  end subroutine allocate_EAMType
390 >
391 >
392 >  subroutine deallocate_EAMType(thisEAMType)
393 >    type (EAMtype), pointer :: thisEAMType
394 >
395 >    ! free Arrays in reverse order of allocation...
396 >    deallocate(thisEAMType%eam_phi_r_pp)      
397 >    deallocate(thisEAMType%eam_rho_r_pp)  
398 >    deallocate(thisEAMType%eam_Z_r_pp)  
399 >    deallocate(thisEAMType%eam_F_rho_pp)  
400 >    deallocate(thisEAMType%eam_phi_r)      
401 >    deallocate(thisEAMType%eam_rho_r)      
402 >    deallocate(thisEAMType%eam_Z_r)  
403 >    deallocate(thisEAMType%eam_F_rho)
404 >    deallocate(thisEAMType%eam_rhovals)
405 >    deallocate(thisEAMType%eam_rvals)
406 >  
407 >  end subroutine deallocate_EAMType
408 >
409 > !! Calculates rho_r
410 >  subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
411 >    integer :: atom1,atom2
412 >    real(kind = dp), dimension(3) :: d
413 >    real(kind = dp), intent(inout)               :: r
414 >    real(kind = dp), intent(inout)               :: rijsq
415 >    ! value of electron density rho do to atom i at atom j
416 >    real(kind = dp) :: rho_i_at_j
417 >    ! value of electron density rho do to atom j at atom i
418 >    real(kind = dp) :: rho_j_at_i
419 >
420 >    ! we don't use the derivatives, dummy variables
421 >    real( kind = dp) :: drho,d2rho
422 >    integer :: eam_err
423 >  
424 >    if (cleanme) call clean_EAM
425 >    cleanme = .false.
426 >
427 >    call  calc_eam_rho(r,rho_i_at_j,drho,d2rho,atom1)
428 >
429 > #ifdef  MPI
430 >    rho_col(atom2) = rho_col(atom2) + rho_i_at_j
431 > #else
432 >    rho(atom2) = rho(atom2) + rho_i_at_j
433 > #endif
434 >
435 >    call calc_eam_rho(r,rho_j_at_i,drho,d2rho,atom2)
436 >
437 > #ifdef  MPI
438 >    rho_row(atom1) = rho_row(atom1) + rho_j_at_i
439 > #else
440 >    rho(atom1) = rho(atom1) + rho_j_at_i
441 > #endif
442 >  end subroutine calc_eam_prepair_rho
443 >
444 >  !! Calculate the functional F(rho) for all atoms
445 >  subroutine calc_eam_prepair_Frho(nlocal,pot)
446 >    integer :: nlocal
447 >    real(kind=dp) :: pot
448 >    integer :: i,j
449 >    real(kind=dp) :: U,U1,U2
450 >    integer :: atype1
451 >    ! reset clean forces to be true at top of calc rho.
452 >    cleanme = .true.
453 >
454 > !! Scatter the electron density in pre-pair
455 > #ifdef MPI
456 >    call scatter(rho_row,rho,plan_row,eam_err)
457 >    if (eam_err /= 0 ) then
458 >      write(errMsg,*) " Error scattering rho_row into rho"
459 >      call handleError(RoutineName,errMesg)
460 >   endif      
461 >    call scatter(rho_col,rho,plan_col,eam_err)
462 >    if (eam_err /= 0 ) then
463 >      write(errMsg,*) " Error scattering rho_col into rho"
464 >      call handleError(RoutineName,errMesg)
465 >   endif
466 > #endif
467 >
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
475 >
476 > #ifdef MPI
477 >    !! communicate f(rho) and derivatives
478 >    call gather(frho,frho_row,plan_row, eam_err)
479 >    if (eam_err /=  0) then
480 >       call handleError("cal_eam_forces()","MPI gather frho_row failure")
481 >    endif
482 >    call gather(dfrhodrho,dfrhodrho_row,plan_row, eam_err)
483 >    if (eam_err /=  0) then
484 >       call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure")
485 >    endif
486 >    call gather(frho,frho_col,plan_col, eam_err)
487 >    if (eam_err /=  0) then
488 >       call handleError("cal_eam_forces()","MPI gather frho_col failure")
489 >    endif
490 >    call gather(dfrhodrho,dfrhodrho_col,plan_col, eam_err)
491 >    if (eam_err /=  0) then
492 >       call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
493 >    endif
494  
495 +    if (nmflag) then
496 +       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_row)
497 +       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_col)
498 +    endif
499 + #endif
500  
501 +
502 +  end subroutine calc_eam_prepair_Frho
503 +
504 +
505 +
506 +
507    subroutine calc_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress)
508   !Arguments    
509      integer, intent(in) ::  atom1, atom2
# Line 258 | Line 513 | contains
513      real( kind = dp ), intent(in), dimension(3) :: d
514      logical, intent(in) :: do_pot, do_stress
515  
516 +    real( kind = dp ) :: drdx,drdy,drdz
517 +    real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
518 +    real( kind = dp ) :: rha,drha,d2rha, dpha
519 +    real( kind = dp ) :: rhb,drhb,d2rhb, dphb
520 +    real( kind = dp ) :: dudr
521 +    real( kind = dp ) :: rci,rcj
522 +    real( kind = dp ) :: drhoidr,drhojdr
523 +    real( kind = dp ) :: d2rhoidrdr
524 +    real( kind = dp ) :: d2rhojdrdr
525 +    real( kind = dp ) :: Fx,Fy,Fz
526 +    real( kind = dp ) :: r,d2pha,phb,d2phb
527 +    integer :: id1,id2
528 +    integer  :: atype1,atype2
529 +
530 +
531   !Local Variables
532      
533  
534 +    
535 +    phab = 0.0E0_DP
536 +    dvpdr = 0.0E0_DP
537 +    d2vpdrdr = 0.0E0_DP
538 +      
539 +
540      if (rij .lt. EAM_rcut) then
541 + #ifdef IS_MPI
542 + !!!!! FIX ME
543 +       atype1 = atid_row(atom1)
544 + #else
545 +       atype1 = atid(atom1)
546 + #endif
547  
548 <       r = dsqrt(rijsq)
549 <       efr(1,j) = -rxij
550 <       efr(2,j) = -ryij
551 <       efr(3,j) = -rzij
552 <
271 <
548 >       drdx = d(1)/rij
549 >       drdy = d(2)/rij
550 >       drdz = d(3)/rij
551 >      
552 >      
553         call calc_eam_rho(r, rha, drha, d2rha, atype1)
554         call calc_eam_phi(r, pha, dpha, d2pha, atype1)
555 <       rci = eam_rcut(eam_atype_map(atype1))
556 < #ifdef MPI
557 <       atype2 = ident_col(j)
555 >  !     rci = eam_rcut(eam_atype_map(atom1))
556 > #ifdef IS_MPI
557 >       atype2 = atid_col(atom2)
558   #else
559 <       atype2 = ident(j)
559 >       atype2 = atid(atom2)
560   #endif
561 <
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))
564 >  !     rcj = eam_rcut(eam_atype_map(atype2))
565  
285 –       phab = 0.0E0_DP
286 –       dvpdr = 0.0E0_DP
287 –       d2vpdrdr = 0.0E0_DP
288 –
566         if (r.lt.rci) then
567            phab = phab + 0.5E0_DP*(rhb/rha)*pha
568            dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
# Line 295 | Line 572 | contains
572                 pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
573                 (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
574         endif
575 +      
576  
299 –
577         if (r.lt.rcj) then
578            phab = phab + 0.5E0_DP*(rha/rhb)*phb
579            dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
# Line 306 | Line 583 | contains
583                 phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
584                 (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
585         endif
586 <
310 <
311 < #ifdef MPI
312 <
313 <       e_row(i) = e_row(i) + phab*0.5
314 <       e_col(i) = e_col(i) + phab*0.5
315 < #else
316 <       if (do_pot) pot = pot + phab
317 < #endif
318 <
586 >      
587         drhoidr = drha
588         drhojdr = drhb
589  
590         d2rhoidrdr = d2rha
591         d2rhojdrdr = d2rhb
592 +
593 +
594   #ifdef MPI
595 <       dudr = drhojdr*dfrhodrho_row(i)+drhoidr*dfrhodrho_col(j) &
595 >       dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
596              + dvpdr
597 <
328 <       if (nmflag) then
329 <          d2 = d2vpdrdr + &
330 <               d2rhoidrdr*dfrhodrho_col(j) + &
331 <               d2rhojdrdr*dfrhodrho_row(i) + &
332 <               drhoidr*drhoidr*d2frhodrhodrho_col(j) + &
333 <               drhojdr*drhojdr*d2frhodrhodrho_row(i)
334 <       endif
597 >    
598   #else
599 <       dudr = drhojdr*dfrhodrho(i)+drhoidr*dfrhodrho(j) &
599 >       dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
600              + dvpdr
601  
339 –       d2 = d2vpdrdr + &
340 –            d2rhoidrdr*dfrhodrho(j) + &
341 –            d2rhojdrdr*dfrhodrho(i) + &
342 –            drhoidr*drhoidr*d2frhodrhodrho(j) + &
343 –            drhojdr*drhojdr*d2frhodrhodrho(i)
602   #endif
603  
604 +       fx = dudr * drdx
605 +       fy = dudr * drdy
606 +       fz = dudr * drdz
607  
347 –       do dim = 1, 3                        
608  
349 –          drdx1 = efr(dim,j) / r
350 –          ftmp = dudr * drdx1
351 –
609   #ifdef MPI
610 <          f_col(dim,j) = f_col(dim,j) - ftmp
611 <          f_row(dim,i) = f_row(dim,i) + ftmp
612 < #else
613 <          f(dim,j) = f(dim,j) - ftmp
357 <          f(dim,i) = f(dim,i) + ftmp
358 < #endif
610 >       if (do_pot) then
611 >          pot_Row(atom1) = pot_Row(atom1) + phab*0.5
612 >          pot_Col(atom2) = pot_Col(atom2) + phab*0.5
613 >       end if
614  
615 <          if (nmflag) then
616 <             idim = 3 * (i-1) + dim
617 <             jdim = 3 * (j-1) + dim
615 >       f_Row(1,atom1) = f_Row(1,atom1) + fx
616 >       f_Row(2,atom1) = f_Row(2,atom1) + fy
617 >       f_Row(3,atom1) = f_Row(3,atom1) + fz
618  
619 <             do dim2 = 1, 3
619 >       f_Col(1,atom2) = f_Col(1,atom2) - fx
620 >       f_Col(2,atom2) = f_Col(2,atom2) - fy
621 >       f_Col(3,atom2) = f_Col(3,atom2) - fz
622 > #else
623 >       if(do_pot) pot = pot + phab
624  
625 <                kt1 = d2 * efr(dim,j) * efr(dim2,j)/r/r
626 <                kt2 = - dudr * efr(dim,j) * efr(dim2,j)/r/r/r
625 >       f(1,atom1) = f(1,atom1) + fx
626 >       f(2,atom1) = f(2,atom1) + fy
627 >       f(3,atom1) = f(3,atom1) + fz
628  
629 <                if (dim.eq.dim2) then
630 <                   kt3 = dudr / r
631 <                else
632 <                   kt3 = 0.0E0_DP
373 <                endif
629 >       f(1,atom2) = f(1,atom2) - fx
630 >       f(2,atom2) = f(2,atom2) - fy
631 >       f(3,atom2) = f(3,atom2) - fz
632 > #endif
633  
375 –                ! The factor of 2 below is to compensate for
376 –                ! overcounting.
377 –                ! Mass weighting is done separately...
634  
379 –                ktmp = (kt1+kt2+kt3)/2.0E0_DP
380 –                idim2 = 3 * (i-1) + dim2
381 –                jdim2 = 3 * (j-1) + dim2
635  
383 –                d(idim,  idim2) = d(idim,idim2)  + ktmp
384 –                d(idim2, idim) = d(idim2,idim)   + ktmp
636  
637 <                d(idim,  jdim2) = d(idim,jdim2)  - ktmp
387 <                d(idim2, jdim) = d(idim2,jdim)   - ktmp
637 >       if (do_stress) then
638  
639 <                d(jdim,  idim2) = d(jdim,idim2)  - ktmp
640 <                d(jdim2, idim) = d(jdim2,idim)   - ktmp
639 > #ifdef MPI
640 >          id1 = tagRow(atom1)
641 >          id2 = tagColumn(atom2)
642 > #else
643 >          id1 = atom1
644 >          id2 = atom2
645 > #endif
646  
647 <                d(jdim,  jdim2) = d(jdim,jdim2)  + ktmp
648 <                d(jdim2, jdim) = d(jdim2,jdim)   + ktmp
647 >          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
648 >            
649 >             tau_Temp(1) = tau_Temp(1) + fx * d(1)
650 >             tau_Temp(2) = tau_Temp(2) + fx * d(2)
651 >             tau_Temp(3) = tau_Temp(3) + fx * d(3)
652 >             tau_Temp(4) = tau_Temp(4) + fy * d(1)
653 >             tau_Temp(5) = tau_Temp(5) + fy * d(2)
654 >             tau_Temp(6) = tau_Temp(6) + fy * d(3)
655 >             tau_Temp(7) = tau_Temp(7) + fz * d(1)
656 >             tau_Temp(8) = tau_Temp(8) + fz * d(2)
657 >             tau_Temp(9) = tau_Temp(9) + fz * d(3)
658 >             virial_Temp = virial_Temp + &
659 >                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
660  
395 –             enddo
661            endif
662 <       enddo
663 <
662 >       endif
663 >      
664      endif
400 – enddo
401 – endif
665  
403 – enddo
666  
667 + end subroutine calc_eam_pair
668  
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  
740 < end subroutine calc_eam_pair
740 >  subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
741  
410 – subroutine calc_eam_rho(r, rho, drho, d2rho, atype)
742  
743 <  !  include 'headers/sizes.h'
743 >    integer :: atype, nx, j
744 >    real( kind = DP ), dimension(:) :: xa
745 >    real( kind = DP ), dimension(:) :: ya
746 >    real( kind = DP ), dimension(:) :: yppa
747 >    real( kind = DP ) :: x, y, dy, d2y
748 >    real( kind = DP ) :: del, h, a, b, c, d
749  
750  
415 – integer atype, etype, number_r
416 – real( kind = DP )  :: r, rho, drho, d2rho
417 – integer :: i
751  
752  
753 < etype = eam_atype_map(atype)
753 >    
754 >    ! this spline code assumes that the x points are equally spaced
755 >    ! do not attempt to use this code if they are not.
756 >    
757 >    
758 >    ! find the closest point with a value below our own:
759 >    j = FLOOR(dble(nx-1) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
760  
761 < if (r.lt.eam_rcut(etype)) then
762 < number_r = eam_nr(etype)
763 < call eam_splint(etype, number_r, eam_rvals, eam_rho_r, &
764 <   eam_rho_r_pp, r, rho, drho, d2rho)
765 < else
766 < rho = 0.0E0_DP
767 < drho = 0.0E0_DP
768 < d2rho = 0.0E0_DP
769 < endif
761 >    ! check to make sure we're inside the spline range:
762 >    if ((j.gt.nx).or.(j.lt.1)) then
763 >       write(default_error,*) "EAM_splint: x is outside bounds of spline"
764 >    endif
765 >    ! check to make sure we haven't screwed up the calculation of j:
766 >    if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
767 >       if (j.ne.nx) then
768 >        write(default_error,*) "EAM_splint: x is outside bounding range"
769 >       endif
770 >    endif
771  
772 < return
773 < end subroutine calc_eam_rho
772 >    del = xa(j+1) - x
773 >    h = xa(j+1) - xa(j)
774 >    
775 >    a = del / h
776 >    b = 1.0E0_DP - a
777 >    c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
778 >    d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
779 >    
780 >    y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
781 >    
782 >    dy = (ya(j+1)-ya(j))/h &
783 >         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
784 >         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
785 >    
786 >    d2y = a*yppa(j) + b*yppa(j+1)
787  
788 < subroutine calc_eam_frho(dens, u, u1, u2, atype)
788 >  end subroutine eam_splint
789  
790 <  ! include 'headers/sizes.h'
790 >  subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
791  
792 < integer atype, etype, number_rho
440 < real( kind = DP ) :: dens, u, u1, u2
441 < real( kind = DP ) :: rho_vals
792 >  
793  
443 – etype = eam_atype_map(atype)
444 – number_rho = eam_nrho(etype)
445 – if (dens.lt.eam_rhovals(number_rho, etype)) then
446 – call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
447 –   eam_f_rho_pp, dens, u, u1, u2)
448 – else
449 – rho_vals = eam_rhovals(number_rho,etype)
450 – call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
451 –   eam_f_rho_pp, rho_vals, u, u1, u2)
452 – endif
794  
795 < return
796 < end subroutine calc_eam_frho
795 >    ! yp1 and ypn are the first derivatives of y at the two endpoints
796 >    ! if boundary is 'L' the lower derivative is used
797 >    ! if boundary is 'U' the upper derivative is used
798 >    ! if boundary is 'B' then both derivatives are used
799 >    ! if boundary is anything else, then both derivatives are assumed to be 0
800 >    
801 >    integer :: nx, i, k, max_array_size
802 >    
803 >    real( kind = DP ), dimension(:)        :: xa
804 >    real( kind = DP ), dimension(:)        :: ya
805 >    real( kind = DP ), dimension(:)        :: yppa
806 >    real( kind = DP ), dimension(size(xa)) :: u
807 >    real( kind = DP ) :: yp1,ypn,un,qn,sig,p
808 >    character boundary
809 >    
810 >    
811 >    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
812 >         (boundary.eq.'b').or.(boundary.eq.'B')) then
813 >       yppa(1) = -0.5E0_DP
814 >       u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-&
815 >            ya(1))/(xa(2)-xa(1))-yp1)
816 >    else
817 >       yppa(1) = 0.0E0_DP
818 >       u(1)  = 0.0E0_DP
819 >    endif
820 >    
821 >    do i = 2, nx - 1
822 >       sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
823 >       p = sig * yppa(i-1) + 2.0E0_DP
824 >       yppa(i) = (sig - 1.0E0_DP) / p
825 >       u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - &
826 >            (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
827 >            (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
828 >    enddo
829 >    
830 >    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
831 >         (boundary.eq.'b').or.(boundary.eq.'B')) then
832 >       qn = 0.5E0_DP
833 >       un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* &
834 >            (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1)))
835 >    else
836 >       qn = 0.0E0_DP
837 >       un = 0.0E0_DP
838 >    endif
839  
840 < subroutine calc_eam_phi(r, phi, dphi, d2phi, atype)
840 >    yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
841 >    
842 >    do k = nx-1, 1, -1
843 >       yppa(k)=yppa(k)*yppa(k+1)+u(k)
844 >    enddo
845  
846 +  end subroutine eam_spline
847  
848  
849  
462 – integer atype, etype, number_r
463 – real( kind = DP ) :: r, phi, dphi, d2phi
850  
851 < etype = eam_atype_map(atype)
466 <
467 < if (r.lt.eam_rcut(etype)) then
468 < number_r = eam_nr(etype)
469 < call eam_splint(etype, number_r, eam_rvals, eam_phi_r, &
470 <   eam_phi_r_pp, r, phi, dphi, d2phi)
471 < else
472 < phi = 0.0E0_DP
473 < dphi = 0.0E0_DP
474 < d2phi = 0.0E0_DP
475 < endif
476 <
477 < return
478 < end subroutine calc_eam_phi
479 <
480 <
481 < subroutine eam_splint(atype, nx, xa, ya, yppa, x, y, dy, d2y)
482 <
483 <  !  include 'headers/sizes.h'
484 <
485 < real( kind = DP ), dimension(:,:) :: xa
486 < real( kind = DP ), dimension(:,:) :: ya
487 < real( kind = DP ), dimension(:,:) :: yppa
488 < real( kind = DP ) :: x, y, dy, d2y
489 < real( kind = DP ) :: del, h, a, b, c, d
490 <
491 <
492 < integer atype, nx, j
493 <
494 <
495 < ! this spline code assumes that the x points are equally spaced
496 < ! do not attempt to use this code if they are not.
497 <
498 <
499 < ! find the closest point with a value below our own:
500 < j = FLOOR(dble(nx-1) * (x - xa(1,atype)) / (xa(nx,atype) - xa(1,atype))) + 1
501 <
502 < ! check to make sure we're inside the spline range:
503 < if ((j.gt.nx).or.(j.lt.1)) call error('eam_splint', &
504 < 'x is outside bounds of spline')
505 <
506 < ! check to make sure we haven't screwed up the calculation of j:
507 < if ((x.lt.xa(j,atype)).or.(x.gt.xa(j+1,atype))) then
508 < if (j.ne.nx) then
509 < call error('eam_splint', &
510 <      'x is outside bounding range')
511 < endif
512 < endif
513 <
514 < del = xa(j+1,atype) - x
515 < h = xa(j+1,atype) - xa(j,atype)
516 <
517 < a = del / h
518 < b = 1.0E0_DP - a
519 < c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
520 < d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
521 <
522 < y = a*ya(j,atype) + b*ya(j+1,atype) + c*yppa(j,atype) + d*yppa(j+1,atype)
523 <
524 < dy = (ya(j+1,atype)-ya(j,atype))/h &
525 < - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j,atype)/6.0E0_DP &
526 < + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1,atype)/6.0E0_DP
527 <
528 < d2y = a*yppa(j,atype) + b*yppa(j+1,atype)
529 <
530 < return
531 < end subroutine eam_splint
532 <
533 < subroutine eam_spline(atype, nx, xa, ya, yppa, yp1, ypn, boundary)
534 <
535 <  !  include 'headers/sizes.h'
536 <
537 <
538 <  ! yp1 and ypn are the first derivatives of y at the two endpoints
539 <  ! if boundary is 'L' the lower derivative is used
540 <  ! if boundary is 'U' the upper derivative is used
541 <  ! if boundary is 'B' then both derivatives are used
542 <  ! if boundary is anything else, then both derivatives are assumed to be 0
543 <
544 < integer nx, i, k, atype, max_array_size
545 <
546 < real( kind = DP ), dimension(:,:) :: xa
547 < real( kind = DP ), dimension(:,:) :: ya
548 < real( kind = DP ), dimension(:,:) :: yppa
549 < real( kind = DP ), allocatable, dimension(:) :: u
550 < real( kind = DP ) :: yp1,ypn,un,qn,sig,p
551 < character boundary
552 <
553 < max_array_size = size(xa,1)
554 < allocate(u(max_array_size))
555 <
556 <
557 < if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
558 < (boundary.eq.'b').or.(boundary.eq.'B')) then
559 < yppa(1, atype) = -0.5E0_DP
560 < u(1) = (3.0E0_DP/(xa(2,atype)-xa(1,atype)))*((ya(2,atype)-&
561 <   ya(1,atype))/(xa(2,atype)-xa(1,atype))-yp1)
562 < else
563 < yppa(1,atype) = 0.0E0_DP
564 < u(1)  = 0.0E0_DP
565 < endif
566 <
567 < do i = 2, nx - 1
568 < sig = (xa(i,atype) - xa(i-1,atype)) / (xa(i+1,atype) - xa(i-1,atype))
569 < p = sig * yppa(i-1,atype) + 2.0E0_DP
570 < yppa(i,atype) = (sig - 1.0E0_DP) / p
571 < u(i) = (6.0E0_DP*((ya(i+1,atype)-ya(i,atype))/(xa(i+1,atype)-xa(i,atype)) - &
572 <   (ya(i,atype)-ya(i-1,atype))/(xa(i,atype)-xa(i-1,atype)))/ &
573 <   (xa(i+1,atype)-xa(i-1,atype)) - sig * u(i-1))/p
574 < enddo
575 <
576 < if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
577 < (boundary.eq.'b').or.(boundary.eq.'B')) then
578 < qn = 0.5E0_DP
579 < un = (3.0E0_DP/(xa(nx,atype)-xa(nx-1,atype)))* &
580 <   (ypn-(ya(nx,atype)-ya(nx-1,atype))/(xa(nx,atype)-xa(nx-1,atype)))
581 < else
582 < qn = 0.0E0_DP
583 < un = 0.0E0_DP
584 < endif
585 <
586 < yppa(nx,atype)=(un-qn*u(nx-1))/(qn*yppa(nx-1,atype)+1.0E0_DP)
587 <
588 < do k = nx-1, 1, -1
589 < yppa(k,atype)=yppa(k,atype)*yppa(k+1,atype)+u(k)
590 < enddo
591 <
592 < deallocate(u)
593 < return
594 < end subroutine eam_spline
595 <
596 <
597 <
598 <
599 < end module calc_eam
851 > end module eam

Diff Legend

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