| 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 | 
| 18 | > | logical, save :: haveRcut = .false. | 
| 19 |  |  | 
| 20 |  | character(len = statusMsgSize) :: errMesg | 
| 21 |  | integer :: eam_err | 
| 51 |  |  | 
| 52 |  |  | 
| 53 |  | !! Arrays for derivatives used in force calculation | 
| 54 | < | real( kind = dp), dimension(:), allocatable :: frho | 
| 55 | < | real( kind = dp), dimension(:), allocatable :: rho | 
| 54 | > | real( kind = dp),save, dimension(:), allocatable :: frho | 
| 55 | > | real( kind = dp),save, dimension(:), allocatable :: rho | 
| 56 |  |  | 
| 57 | < | real( kind = dp), dimension(:), allocatable :: dfrhodrho | 
| 58 | < | real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho | 
| 57 | > | real( kind = dp),save, dimension(:), allocatable :: dfrhodrho | 
| 58 | > | real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho | 
| 59 |  |  | 
| 60 |  |  | 
| 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 | 
| 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 :: d2frhodrhodrho_col | 
| 70 | > | real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row | 
| 71 |  | #endif | 
| 72 |  |  | 
| 73 |  | type, private :: EAMTypeList | 
| 117 |  |  | 
| 118 |  | status = 0 | 
| 119 |  |  | 
| 120 | < | write(*,*) "Adding new eamtype: ",eam_ident | 
| 120 | > |  | 
| 121 |  | !! Assume that atypes has already been set and get the total number of types in atypes | 
| 122 |  | !! Also assume that every member of atypes is a EAM model. | 
| 123 |  |  | 
| 222 |  | enddo | 
| 223 |  |  | 
| 224 |  |  | 
| 225 | < | current_rcut_max = EAMList%EAMParams(1)%eam_rcut | 
| 226 | < | !! find the smallest rcut for any eam atype | 
| 227 | < | do i = 2, EAMList%currentAddition | 
| 228 | < | current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut) | 
| 229 | < | end do | 
| 225 | > | !       current_rcut_max = EAMList%EAMParams(1)%eam_rcut | 
| 226 | > | !! find the smallest rcut for any eam atype | 
| 227 | > | !       do i = 2, EAMList%currentAddition | 
| 228 | > | !          current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut) | 
| 229 | > | !       end do | 
| 230 |  |  | 
| 231 | < | EAM_rcut = current_rcut_max | 
| 232 | < | EAM_rcut_orig = current_rcut_max | 
| 231 | > | !       EAM_rcut = current_rcut_max | 
| 232 | > | !       EAM_rcut_orig = current_rcut_max | 
| 233 |  | !       do i = 1, EAMList%currentAddition | 
| 234 |  | !          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i | 
| 235 |  | !       end do | 
| 350 |  |  | 
| 351 |  | end subroutine allocateEAM | 
| 352 |  |  | 
| 353 | + | !! C sets rcut to be the largest cutoff of any atype | 
| 354 | + | !! present in this simulation. Doesn't include all atypes | 
| 355 | + | !! sim knows about, just those in the simulation. | 
| 356 |  | subroutine setCutoffEAM(rcut, status) | 
| 357 |  | real(kind=dp) :: rcut | 
| 358 |  | integer :: status | 
| 359 | + | status = 0 | 
| 360 |  |  | 
| 361 | < | if (rcut < EAM_rcut) then | 
| 358 | < | EAM_rcut = rcut | 
| 359 | < | endif | 
| 361 | > | EAM_rcut = rcut | 
| 362 |  |  | 
| 361 | – |  | 
| 363 |  | end subroutine setCutoffEAM | 
| 364 |  |  | 
| 365 |  |  | 
| 366 |  |  | 
| 367 |  | subroutine clean_EAM() | 
| 368 | < |  | 
| 368 | < | ! clean non-IS_MPI first | 
| 368 | > | ! clean non-IS_MPI first | 
| 369 |  | frho = 0.0_dp | 
| 370 |  | rho  = 0.0_dp | 
| 371 |  | dfrhodrho = 0.0_dp | 
| 484 |  | integer :: myid_atom2 | 
| 485 |  |  | 
| 486 |  | ! check to see if we need to be cleaned at the start of a force loop | 
| 487 | – | if (cleanme) call clean_EAM | 
| 488 | – | cleanme = .false. | 
| 487 |  |  | 
| 488 | + | if (cleanme) then | 
| 489 | + | call clean_EAM | 
| 490 | + | cleanme = .false. | 
| 491 | + | end if | 
| 492 | + |  | 
| 493 |  |  | 
| 494 | + |  | 
| 495 | + |  | 
| 496 |  | #ifdef IS_MPI | 
| 497 |  | myid_atom1 = atid_Row(atom1) | 
| 498 |  | myid_atom2 = atid_Col(atom2) | 
| 537 |  | #endif | 
| 538 |  | endif | 
| 539 |  |  | 
| 540 | + |  | 
| 541 | + |  | 
| 542 |  | end subroutine calc_eam_prepair_rho | 
| 543 |  |  | 
| 544 |  |  | 
| 556 |  | integer :: n_rho_points | 
| 557 |  | ! reset clean forces to be true at top of calc rho. | 
| 558 |  | cleanme = .true. | 
| 559 | < |  | 
| 559 | > |  | 
| 560 |  | !! Scatter the electron density from  pre-pair calculation back to local atoms | 
| 561 |  | #ifdef IS_MPI | 
| 562 |  | call scatter(rho_row,rho,plan_row,eam_err) | 
| 601 |  | pot = pot + u | 
| 602 |  | enddo | 
| 603 |  |  | 
| 604 | + |  | 
| 605 |  |  | 
| 598 | – |  | 
| 606 |  | #ifdef IS_MPI | 
| 607 |  | !! communicate f(rho) and derivatives back into row and column arrays | 
| 608 |  | call gather(frho,frho_row,plan_row, eam_err) | 
| 668 |  | !Local Variables | 
| 669 |  |  | 
| 670 |  |  | 
| 671 | < |  | 
| 671 | > |  | 
| 672 |  | phab = 0.0E0_DP | 
| 673 |  | dvpdr = 0.0E0_DP | 
| 674 |  | d2vpdrdr = 0.0E0_DP | 
| 675 | < |  | 
| 669 | < |  | 
| 675 | > |  | 
| 676 |  | if (rij .lt. EAM_rcut) then | 
| 677 |  | #ifdef IS_MPI | 
| 678 |  | !!!!! FIX ME | 
| 759 |  | #ifdef IS_MPI | 
| 760 |  | dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) & | 
| 761 |  | + dvpdr | 
| 762 | < |  | 
| 762 | > |  | 
| 763 |  | #else | 
| 764 |  | dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) & | 
| 765 |  | + dvpdr | 
| 766 | < |  | 
| 766 | > | ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2) | 
| 767 |  | #endif | 
| 768 |  |  | 
| 769 |  | fx = dudr * drdx | 
| 785 |  | f_Col(2,atom2) = f_Col(2,atom2) - fy | 
| 786 |  | f_Col(3,atom2) = f_Col(3,atom2) - fz | 
| 787 |  | #else | 
| 788 | < | if(do_pot) pot = pot + phab | 
| 789 | < |  | 
| 788 | > |  | 
| 789 | > | if(do_pot) then | 
| 790 | > | pot = pot + phab | 
| 791 | > | end if | 
| 792 | > |  | 
| 793 |  | f(1,atom1) = f(1,atom1) + fx | 
| 794 |  | f(2,atom1) = f(2,atom1) + fy | 
| 795 |  | f(3,atom1) = f(3,atom1) + fz |