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

Diff Legend

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