1 |
< |
!! |
1 |
> |
!! |
2 |
|
!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
|
!! |
4 |
|
!! The University of Notre Dame grants you ("Licensee") a |
66 |
|
logical, save :: haveRcut = .false. |
67 |
|
logical, save :: haveMixingMap = .false. |
68 |
|
logical, save :: useGeometricDistanceMixing = .false. |
69 |
< |
|
70 |
< |
|
69 |
> |
logical, save :: cleanArrays = .true. |
70 |
> |
logical, save :: arraysAllocated = .false. |
71 |
|
|
72 |
|
|
73 |
|
character(len = statusMsgSize) :: errMesg |
76 |
|
character(len = 200) :: errMsg |
77 |
|
character(len=*), parameter :: RoutineName = "Sutton-Chen MODULE" |
78 |
|
!! Logical that determines if eam arrays should be zeroed |
79 |
– |
logical :: cleanme = .true. |
79 |
|
logical :: nmflag = .false. |
80 |
|
|
81 |
|
|
272 |
|
MixingMap(i,j)%rcut = 2.0_dp *MixingMap(i,j)%alpha |
273 |
|
MixingMap(i,j)%vpair_pot = MixingMap(i,j)%epsilon* & |
274 |
|
(MixingMap(i,j)%alpha/MixingMap(i,j)%rcut)**MixingMap(i,j)%n |
275 |
+ |
|
276 |
|
if (i.ne.j) then |
277 |
|
MixingMap(j,i)%epsilon = MixingMap(i,j)%epsilon |
278 |
|
MixingMap(j,i)%m = MixingMap(i,j)%m |
291 |
|
|
292 |
|
!! routine checks to see if array is allocated, deallocates array if allocated |
293 |
|
!! and then creates the array to the required size |
294 |
< |
subroutine allocateSC(status) |
295 |
< |
integer, intent(out) :: status |
294 |
> |
subroutine allocateSC() |
295 |
> |
integer :: status |
296 |
|
|
297 |
|
#ifdef IS_MPI |
298 |
|
integer :: nAtomsInRow |
300 |
|
#endif |
301 |
|
integer :: alloc_stat |
302 |
|
|
303 |
< |
|
303 |
> |
|
304 |
|
status = 0 |
305 |
|
#ifdef IS_MPI |
306 |
|
nAtomsInRow = getNatomsInRow(plan_atom_row) |
313 |
|
allocate(frho(nlocal),stat=alloc_stat) |
314 |
|
if (alloc_stat /= 0) then |
315 |
|
status = -1 |
316 |
– |
return |
316 |
|
end if |
317 |
|
|
318 |
|
if (allocated(rho)) deallocate(rho) |
319 |
|
allocate(rho(nlocal),stat=alloc_stat) |
320 |
|
if (alloc_stat /= 0) then |
321 |
|
status = -1 |
323 |
– |
return |
322 |
|
end if |
323 |
|
|
324 |
|
if (allocated(dfrhodrho)) deallocate(dfrhodrho) |
325 |
|
allocate(dfrhodrho(nlocal),stat=alloc_stat) |
326 |
|
if (alloc_stat /= 0) then |
327 |
|
status = -1 |
330 |
– |
return |
328 |
|
end if |
329 |
|
|
330 |
|
#ifdef IS_MPI |
333 |
|
allocate(rho_tmp(nlocal),stat=alloc_stat) |
334 |
|
if (alloc_stat /= 0) then |
335 |
|
status = -1 |
339 |
– |
return |
336 |
|
end if |
337 |
|
|
338 |
|
|
340 |
|
allocate(frho_row(nAtomsInRow),stat=alloc_stat) |
341 |
|
if (alloc_stat /= 0) then |
342 |
|
status = -1 |
347 |
– |
return |
343 |
|
end if |
344 |
|
if (allocated(rho_row)) deallocate(rho_row) |
345 |
|
allocate(rho_row(nAtomsInRow),stat=alloc_stat) |
346 |
|
if (alloc_stat /= 0) then |
347 |
|
status = -1 |
353 |
– |
return |
348 |
|
end if |
349 |
|
if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row) |
350 |
|
allocate(dfrhodrho_row(nAtomsInRow),stat=alloc_stat) |
351 |
|
if (alloc_stat /= 0) then |
352 |
|
status = -1 |
359 |
– |
return |
353 |
|
end if |
354 |
|
|
355 |
|
|
359 |
|
allocate(frho_col(nAtomsInCol),stat=alloc_stat) |
360 |
|
if (alloc_stat /= 0) then |
361 |
|
status = -1 |
369 |
– |
return |
362 |
|
end if |
363 |
|
if (allocated(rho_col)) deallocate(rho_col) |
364 |
|
allocate(rho_col(nAtomsInCol),stat=alloc_stat) |
365 |
|
if (alloc_stat /= 0) then |
366 |
|
status = -1 |
375 |
– |
return |
367 |
|
end if |
368 |
|
if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col) |
369 |
|
allocate(dfrhodrho_col(nAtomsInCol),stat=alloc_stat) |
370 |
|
if (alloc_stat /= 0) then |
371 |
|
status = -1 |
381 |
– |
return |
372 |
|
end if |
373 |
|
|
374 |
|
#endif |
375 |
< |
|
375 |
> |
if (status == -1) then |
376 |
> |
call handleError("SuttonChen:allocateSC","Error in allocating SC arrays") |
377 |
> |
end if |
378 |
> |
arraysAllocated = .true. |
379 |
|
end subroutine allocateSC |
380 |
|
|
381 |
|
!! C sets rcut to be the largest cutoff of any atype |
390 |
|
|
391 |
|
end subroutine setCutoffSC |
392 |
|
|
393 |
+ |
!! This array allocates module arrays if needed and builds mixing map. |
394 |
|
subroutine clean_SC() |
395 |
< |
|
395 |
> |
if (.not.arraysAllocated) call allocateSC() |
396 |
> |
if (.not.haveMixingMap) call createMixingMap() |
397 |
|
! clean non-IS_MPI first |
398 |
|
frho = 0.0_dp |
399 |
|
rho = 0.0_dp |
434 |
|
|
435 |
|
! check to see if we need to be cleaned at the start of a force loop |
436 |
|
|
437 |
+ |
if (cleanArrays) call clean_SC() |
438 |
+ |
cleanArrays = .false. |
439 |
|
|
443 |
– |
|
444 |
– |
|
440 |
|
#ifdef IS_MPI |
441 |
|
Atid1 = Atid_row(Atom1) |
442 |
|
Atid2 = Atid_col(Atom2) |
471 |
|
real(kind=dp) :: pot |
472 |
|
integer :: i,j |
473 |
|
integer :: atom |
479 |
– |
real(kind=dp) :: U,U1,U2 |
474 |
|
integer :: atype1 |
475 |
|
integer :: atid1 |
476 |
|
integer :: myid |
477 |
|
|
478 |
|
|
485 |
– |
cleanme = .true. |
479 |
|
!! Scatter the electron density from pre-pair calculation back to local atoms |
480 |
|
#ifdef IS_MPI |
481 |
|
call scatter(rho_row,rho,plan_atom_row,sc_err) |
499 |
|
do atom = 1, nlocal |
500 |
|
Myid = SCList%atidtoSctype(Atid(atom)) |
501 |
|
frho(atom) = -SCList%SCTypes(Myid)%c * & |
502 |
< |
SCList%SCTypes(Myid)%epsilon * sqrt(rho(i)) |
502 |
> |
SCList%SCTypes(Myid)%epsilon * sqrt(rho(atom)) |
503 |
|
|
504 |
|
dfrhodrho(atom) = 0.5_dp*frho(atom)/rho(atom) |
505 |
< |
pot = pot + u |
505 |
> |
|
506 |
> |
pot = pot + frho(atom) |
507 |
|
enddo |
508 |
|
|
509 |
|
#ifdef IS_MPI |
547 |
|
real( kind = dp ) :: dvpdr |
548 |
|
real( kind = dp ) :: drhodr |
549 |
|
real( kind = dp ) :: dudr |
556 |
– |
real( kind = dp ) :: rcij |
550 |
|
real( kind = dp ) :: drhoidr,drhojdr |
551 |
|
real( kind = dp ) :: Fx,Fy,Fz |
552 |
< |
real( kind = dp ) :: r,d2pha,phb,d2phb |
552 |
> |
real( kind = dp ) :: d2pha,phb,d2phb |
553 |
|
real( kind = dp ) :: pot_temp,vptmp |
554 |
|
real( kind = dp ) :: epsilonij,aij,nij,mij,vcij |
555 |
|
integer :: id1,id2 |
559 |
|
!Local Variables |
560 |
|
|
561 |
|
! write(*,*) "Frho: ", Frho(atom1) |
562 |
+ |
|
563 |
+ |
cleanArrays = .true. |
564 |
|
|
570 |
– |
|
565 |
|
dvpdr = 0.0E0_DP |
566 |
|
|
567 |
|
|
587 |
|
mij = MixingMap(mytype_atom1,mytype_atom2)%m |
588 |
|
vcij = MixingMap(mytype_atom1,mytype_atom2)%vpair_pot |
589 |
|
|
590 |
< |
vptmp = epsilonij*((aij/r)**nij) |
590 |
> |
vptmp = epsilonij*((aij/rij)**nij) |
591 |
|
|
592 |
|
|
593 |
< |
dvpdr = -nij*vptmp/r |
594 |
< |
drhodr = -mij*((aij/r)**mij)/r |
593 |
> |
dvpdr = -nij*vptmp/rij |
594 |
> |
drhodr = -mij*((aij/rij)**mij)/rij |
595 |
|
|
596 |
|
|
597 |
|
dudr = drhodr*(dfrhodrho(atom1)+dfrhodrho(atom2)) & |