ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE/libmdtools/calc_eam.F90

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

Diff Legend

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