ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2432 by chuckv, Tue Nov 15 16:01:06 2005 UTC vs.
Revision 2717 by gezelter, Mon Apr 17 21:49:12 2006 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.68 2005-11-15 16:01:06 chuckv Exp $, $Date: 2005-11-15 16:01:06 $, $Name: not supported by cvs2svn $, $Revision: 1.68 $
48 > !! @version $Id: doForces.F90,v 1.78 2006-04-17 21:49:12 gezelter Exp $, $Date: 2006-04-17 21:49:12 $, $Name: not supported by cvs2svn $, $Revision: 1.78 $
49  
50  
51   module doForces
# Line 64 | Line 64 | module doForces
64    use eam
65    use suttonchen
66    use status
67 +  use interpolation
68   #ifdef IS_MPI
69    use mpiSimulation
70   #endif
# Line 87 | Line 88 | module doForces
88    logical, save :: haveInteractionHash = .false.
89    logical, save :: haveGtypeCutoffMap = .false.
90    logical, save :: haveDefaultCutoffs = .false.
91 <  logical, save :: haveRlist = .false.
91 >  logical, save :: haveSkinThickness = .false.
92 >  logical, save :: haveElectrostaticSummationMethod = .false.
93 >  logical, save :: haveCutoffPolicy = .false.
94 >  logical, save :: VisitCutoffsAfterComputing = .false.
95  
96    logical, save :: FF_uses_DirectionalAtoms
97    logical, save :: FF_uses_Dipoles
# Line 106 | Line 110 | module doForces
110    logical, save :: SIM_uses_PBC
111  
112    integer, save :: electrostaticSummationMethod
113 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
114  
115 +  real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
116 +  real(kind=dp), save :: skinThickness
117 +  logical, save :: defaultDoShift
118 +
119    public :: init_FF
120 <  public :: setDefaultCutoffs
120 >  public :: setCutoffs
121 >  public :: cWasLame
122 >  public :: setElectrostaticMethod
123 >  public :: setCutoffPolicy
124 >  public :: setSkinThickness
125    public :: do_force_loop
113 –  public :: createInteractionHash
114 –  public :: createGtypeCutoffMap
115 –  public :: getStickyCut
116 –  public :: getStickyPowerCut
117 –  public :: getGayBerneCut
118 –  public :: getEAMCut
119 –  public :: getShapeCut
126  
127   #ifdef PROFILE
128    public :: getforcetime
# Line 144 | Line 150 | module doForces
150    end type gtypeCutoffs
151    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
152  
147 –  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
148 –  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
149 –  real(kind=dp),save :: listSkin
150 –  
153   contains
154  
155 <  subroutine createInteractionHash(status)
155 >  subroutine createInteractionHash()
156      integer :: nAtypes
155 –    integer, intent(out) :: status
157      integer :: i
158      integer :: j
159      integer :: iHash
# Line 177 | Line 178 | contains
178      logical :: j_is_MEAM
179      real(kind=dp) :: myRcut
180  
180 –
181 –    status = 0  
182 –
181      if (.not. associated(atypes)) then
182 <       call handleError("atype", "atypes was not present before call of createInteractionHash!")
185 <       status = -1
182 >       call handleError("doForces", "atypes was not present before call of createInteractionHash!")
183         return
184      endif
185      
186      nAtypes = getSize(atypes)
187      
188      if (nAtypes == 0) then
189 <       status = -1
189 >       call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
190         return
191      end if
192  
# Line 276 | Line 273 | contains
273      haveInteractionHash = .true.
274    end subroutine createInteractionHash
275  
276 <  subroutine createGtypeCutoffMap(stat)
276 >  subroutine createGtypeCutoffMap()
277  
281 –    integer, intent(out), optional :: stat
278      logical :: i_is_LJ
279      logical :: i_is_Elect
280      logical :: i_is_Sticky
# Line 286 | Line 282 | contains
282      logical :: i_is_GB
283      logical :: i_is_EAM
284      logical :: i_is_Shape
285 +    logical :: i_is_SC
286      logical :: GtypeFound
287  
288      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
# Line 293 | Line 290 | contains
290      integer :: nGroupsInRow
291      integer :: nGroupsInCol
292      integer :: nGroupTypesRow,nGroupTypesCol
293 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
293 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
294      real(kind=dp) :: biggestAtypeCutoff
295  
299 –    stat = 0
296      if (.not. haveInteractionHash) then
297 <       call createInteractionHash(myStatus)      
302 <       if (myStatus .ne. 0) then
303 <          write(default_error, *) 'createInteractionHash failed in doForces!'
304 <          stat = -1
305 <          return
306 <       endif
297 >       call createInteractionHash()      
298      endif
299   #ifdef IS_MPI
300      nGroupsInRow = getNgroupsInRow(plan_group_row)
# Line 321 | Line 312 | contains
312            call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
313            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
314            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
315 <          
315 >          call getElementProperty(atypes, i, "is_SC", i_is_SC)
316  
317            if (haveDefaultCutoffs) then
318               atypeMaxCutoff(i) = defaultRcut
# Line 354 | Line 345 | contains
345                  thisRcut = getShapeCut(i)
346                  if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
347               endif
348 +             if (i_is_SC) then
349 +                thisRcut = getSCCut(i)
350 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
351 +             endif
352            endif
353 <          
359 <          
353 >                    
354            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
355               biggestAtypeCutoff = atypeMaxCutoff(i)
356            endif
357  
358         endif
359      enddo
366 –  
367 –
360      
361      istart = 1
362      jstart = 1
# Line 408 | Line 400 | contains
400         allocate(groupToGtypeCol(jend))
401      end if
402  
403 <    if(.not.associated(groupToGtypeCol)) then
404 <       allocate(groupToGtypeCol(jend))
403 >    if(.not.associated(groupMaxCutoffCol)) then
404 >       allocate(groupMaxCutoffCol(jend))
405      else
406 <       deallocate(groupToGtypeCol)
407 <       allocate(groupToGtypeCol(jend))
406 >       deallocate(groupMaxCutoffCol)
407 >       allocate(groupMaxCutoffCol(jend))
408      end if
409      if(.not.associated(gtypeMaxCutoffCol)) then
410         allocate(gtypeMaxCutoffCol(jend))
# Line 435 | Line 427 | contains
427      
428      tol = 1.0d-6
429      nGroupTypesRow = 0
430 <
430 >    nGroupTypesCol = 0
431      do i = istart, iend      
432         n_in_i = groupStartRow(i+1) - groupStartRow(i)
433         groupMaxCutoffRow(i) = 0.0_dp
# Line 450 | Line 442 | contains
442               groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
443            endif          
444         enddo
453 –
445         if (nGroupTypesRow.eq.0) then
446            nGroupTypesRow = nGroupTypesRow + 1
447            gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
# Line 512 | Line 503 | contains
503      gtypeMaxCutoffCol => gtypeMaxCutoffRow
504      groupMaxCutoffCol => groupMaxCutoffRow
505   #endif
515 –
516 –
506  
518 –
519 –
507      !! allocate the gtypeCutoffMap here.
508      allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
509      !! then we do a double loop over all the group TYPES to find the cutoff
510      !! map between groups of two types
511      tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
512  
513 <    do i = 1, nGroupTypesRow
513 >    do i = 1, nGroupTypesRow      
514         do j = 1, nGroupTypesCol
515        
516            select case(cutoffPolicy)
# Line 538 | Line 525 | contains
525               return
526            end select
527            gtypeCutoffMap(i,j)%rcut = thisRcut
528 +          
529 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
530 +
531            gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
542 –          skin = defaultRlist - defaultRcut
543 –          listSkin = skin ! set neighbor list skin thickness
544 –          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
532  
533 +          if (.not.haveSkinThickness) then
534 +             skinThickness = 1.0_dp
535 +          endif
536 +
537 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
538 +
539            ! sanity check
540  
541            if (haveDefaultCutoffs) then
# Line 552 | Line 545 | contains
545            endif
546         enddo
547      enddo
548 +
549      if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
550      if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
551      if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
# Line 565 | Line 559 | contains
559      haveGtypeCutoffMap = .true.
560     end subroutine createGtypeCutoffMap
561  
562 <   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
569 <     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
570 <     integer, intent(in) :: cutPolicy
562 >   subroutine setCutoffs(defRcut, defRsw)
563  
564 +     real(kind=dp),intent(in) :: defRcut, defRsw
565 +     character(len = statusMsgSize) :: errMsg
566 +     integer :: localError
567 +
568       defaultRcut = defRcut
569       defaultRsw = defRsw
570 <     defaultRlist = defRlist
571 <     cutoffPolicy = cutPolicy
570 >    
571 >     defaultDoShift = .false.
572 >     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
573 >        
574 >        write(errMsg, *) &
575 >             'cutoffRadius and switchingRadius are set to the same', newline &
576 >             // tab, 'value.  OOPSE will use shifted ', newline &
577 >             // tab, 'potentials instead of switching functions.'
578 >        
579 >        call handleInfo("setCutoffs", errMsg)
580 >        
581 >        defaultDoShift = .true.
582 >        
583 >     endif
584  
585 +     localError = 0
586 +     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
587 +     call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
588 +     call setCutoffEAM( defaultRcut )
589 +     call setCutoffSC( defaultRcut )
590 +     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
591 +     call setHmatDangerousRcutValue(defaultRcut)
592 +
593       haveDefaultCutoffs = .true.
594 <   end subroutine setDefaultCutoffs
594 >     haveGtypeCutoffMap = .false.
595 >   end subroutine setCutoffs
596  
597 +   subroutine cWasLame()
598 +    
599 +     VisitCutoffsAfterComputing = .true.
600 +     return
601 +    
602 +   end subroutine cWasLame
603 +  
604     subroutine setCutoffPolicy(cutPolicy)
605 <
605 >    
606       integer, intent(in) :: cutPolicy
607 +    
608       cutoffPolicy = cutPolicy
609 <     call createGtypeCutoffMap()
609 >     haveCutoffPolicy = .true.
610 >     haveGtypeCutoffMap = .false.
611 >    
612     end subroutine setCutoffPolicy
613 +  
614 +   subroutine setElectrostaticMethod( thisESM )
615 +
616 +     integer, intent(in) :: thisESM
617 +
618 +     electrostaticSummationMethod = thisESM
619 +     haveElectrostaticSummationMethod = .true.
620      
621 +   end subroutine setElectrostaticMethod
622 +
623 +   subroutine setSkinThickness( thisSkin )
624      
625 <  subroutine setSimVariables()
626 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
627 <    SIM_uses_EAM = SimUsesEAM()
628 <    SIM_uses_SC  = SimUsesSC()
629 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
630 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
631 <    SIM_uses_PBC = SimUsesPBC()
625 >     real(kind=dp), intent(in) :: thisSkin
626 >    
627 >     skinThickness = thisSkin
628 >     haveSkinThickness = .true.    
629 >     haveGtypeCutoffMap = .false.
630 >    
631 >   end subroutine setSkinThickness
632 >      
633 >   subroutine setSimVariables()
634 >     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
635 >     SIM_uses_EAM = SimUsesEAM()
636 >     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
637 >     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
638 >     SIM_uses_PBC = SimUsesPBC()
639 >     SIM_uses_SC = SimUsesSC()
640 >    
641 >     haveSIMvariables = .true.
642 >    
643 >     return
644 >   end subroutine setSimVariables
645  
596 –    haveSIMvariables = .true.
597 –
598 –    return
599 –  end subroutine setSimVariables
600 –
646    subroutine doReadyCheck(error)
647      integer, intent(out) :: error
648  
# Line 606 | Line 651 | contains
651      error = 0
652  
653      if (.not. haveInteractionHash) then      
654 <       myStatus = 0      
610 <       call createInteractionHash(myStatus)      
611 <       if (myStatus .ne. 0) then
612 <          write(default_error, *) 'createInteractionHash failed in doForces!'
613 <          error = -1
614 <          return
615 <       endif
654 >       call createInteractionHash()      
655      endif
656  
657      if (.not. haveGtypeCutoffMap) then        
658 <       myStatus = 0      
620 <       call createGtypeCutoffMap(myStatus)      
621 <       if (myStatus .ne. 0) then
622 <          write(default_error, *) 'createGtypeCutoffMap failed in doForces!'
623 <          error = -1
624 <          return
625 <       endif
658 >       call createGtypeCutoffMap()      
659      endif
660 +
661 +
662 +    if (VisitCutoffsAfterComputing) then
663 +       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
664 +       call setHmatDangerousRcutValue(largestRcut)
665 +       call setLJsplineRmax(largestRcut)
666 +       call setCutoffEAM(largestRcut)
667 +       call setCutoffSC(largestRcut)
668 +       VisitCutoffsAfterComputing = .false.
669 +    endif
670 +
671  
672      if (.not. haveSIMvariables) then
673         call setSimVariables()
# Line 658 | Line 702 | contains
702    end subroutine doReadyCheck
703  
704  
705 <  subroutine init_FF(thisESM, thisStat)
705 >  subroutine init_FF(thisStat)
706  
663 –    integer, intent(in) :: thisESM
707      integer, intent(out) :: thisStat  
708      integer :: my_status, nMatches
709      integer, pointer :: MatchList(:) => null()
# Line 668 | Line 711 | contains
711      !! assume things are copacetic, unless they aren't
712      thisStat = 0
713  
671 –    electrostaticSummationMethod = thisESM
672 –
714      !! init_FF is called *after* all of the atom types have been
715      !! defined in atype_module using the new_atype subroutine.
716      !!
# Line 680 | Line 721 | contains
721      FF_uses_Dipoles = .false.
722      FF_uses_GayBerne = .false.
723      FF_uses_EAM = .false.
724 +    FF_uses_SC = .false.
725  
726      call getMatchingElementList(atypes, "is_Directional", .true., &
727           nMatches, MatchList)
# Line 695 | Line 737 | contains
737  
738      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
739      if (nMatches .gt. 0) FF_uses_EAM = .true.
740 +
741 +    call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
742 +    if (nMatches .gt. 0) FF_uses_SC = .true.
743  
744  
745      haveSaneForceField = .true.
# Line 766 | Line 811 | contains
811      real( kind = DP ) :: rVal
812      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
813      real(kind=dp) :: rfpot, mu_i, virial
814 +    real(kind=dp):: rCut
815      integer :: me_i, me_j, n_in_i, n_in_j
816      logical :: is_dp_i
817      integer :: neighborListSize
# Line 839 | Line 885 | contains
885         ! (but only on the first time through):
886         if (loop .eq. loopStart) then
887   #ifdef IS_MPI
888 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
888 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
889                 update_nlist)
890   #else
891 <          call checkNeighborList(nGroups, q_group, listSkin, &
891 >          call checkNeighborList(nGroups, q_group, skinThickness, &
892                 update_nlist)
893   #endif
894         endif
# Line 922 | Line 968 | contains
968  
969                     list(nlist) = j
970                  endif
971 +
972 +
973                  
974                  if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
975  
976 +                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
977                     if (loop .eq. PAIR_LOOP) then
978                        vij = 0.0d0
979 <                      fij(1:3) = 0.0d0
979 >                      fij(1) = 0.0_dp
980 >                      fij(2) = 0.0_dp
981 >                      fij(3) = 0.0_dp
982                     endif
983                    
984 <                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
985 <                        in_switching_region)
984 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
985 >                        group_switch, in_switching_region)
986                    
987                     n_in_j = groupStartCol(j+1) - groupStartCol(j)
988                    
# Line 946 | Line 997 | contains
997                           if (skipThisPair(atom1, atom2))  cycle inner
998                          
999                           if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1000 <                            d_atm(1:3) = d_grp(1:3)
1000 >                            d_atm(1) = d_grp(1)
1001 >                            d_atm(2) = d_grp(2)
1002 >                            d_atm(3) = d_grp(3)
1003                              ratmsq = rgrpsq
1004                           else
1005   #ifdef IS_MPI
# Line 961 | Line 1014 | contains
1014                           if (loop .eq. PREPAIR_LOOP) then
1015   #ifdef IS_MPI                      
1016                              call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1017 <                                 rgrpsq, d_grp, do_pot, do_stress, &
1017 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1018                                   eFrame, A, f, t, pot_local)
1019   #else
1020                              call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1021 <                                 rgrpsq, d_grp, do_pot, do_stress, &
1021 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1022                                   eFrame, A, f, t, pot)
1023   #endif                                              
1024                           else
1025   #ifdef IS_MPI                      
1026                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1027                                   do_pot, eFrame, A, f, t, pot_local, vpair, &
1028 <                                 fpair, d_grp, rgrp)
1028 >                                 fpair, d_grp, rgrp, rCut)
1029   #else
1030                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1031                                   do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1032 <                                 d_grp, rgrp)
1032 >                                 d_grp, rgrp, rCut)
1033   #endif
1034                              vij = vij + vpair
1035 <                            fij(1:3) = fij(1:3) + fpair(1:3)
1035 >                            fij(1) = fij(1) + fpair(1)
1036 >                            fij(2) = fij(2) + fpair(2)
1037 >                            fij(3) = fij(3) + fpair(3)
1038                           endif
1039                        enddo inner
1040                     enddo
# Line 1184 | Line 1239 | contains
1239    end subroutine do_force_loop
1240  
1241    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1242 <       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1242 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1243  
1244      real( kind = dp ) :: vpair, sw
1245      real( kind = dp ), dimension(LR_POT_TYPES) :: pot
# Line 1201 | Line 1256 | contains
1256      real ( kind = dp ), intent(inout) :: r_grp
1257      real ( kind = dp ), intent(inout) :: d(3)
1258      real ( kind = dp ), intent(inout) :: d_grp(3)
1259 +    real ( kind = dp ), intent(inout) :: rCut
1260      real ( kind = dp ) :: r
1261      integer :: me_i, me_j
1262  
# Line 1221 | Line 1277 | contains
1277      iHash = InteractionHash(me_i, me_j)
1278      
1279      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1280 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1280 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1281              pot(VDW_POT), f, do_pot)
1282      endif
1283      
1284      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1285 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1285 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1286              pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1287      endif
1288      
# Line 1246 | Line 1302 | contains
1302      endif
1303      
1304      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1305 <       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1305 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1306              pot(VDW_POT), A, f, t, do_pot)
1307      endif
1308      
# Line 1266 | Line 1322 | contains
1322      endif
1323  
1324      if ( iand(iHash, SC_PAIR).ne.0 ) then      
1325 <       call do_SC_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1325 >       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1326              pot(METALLIC_POT), f, do_pot)
1327      endif
1328  
# Line 1274 | Line 1330 | contains
1330      
1331    end subroutine do_pair
1332  
1333 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1333 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1334         do_pot, do_stress, eFrame, A, f, t, pot)
1335  
1336      real( kind = dp ) :: sw
# Line 1286 | Line 1342 | contains
1342  
1343      logical, intent(inout) :: do_pot, do_stress
1344      integer, intent(in) :: i, j
1345 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1345 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1346      real ( kind = dp )                :: r, rc
1347      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1348  
# Line 1305 | Line 1361 | contains
1361      iHash = InteractionHash(me_i, me_j)
1362  
1363      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1364 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1364 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1365      endif
1366  
1367      if ( iand(iHash, SC_PAIR).ne.0 ) then      
1368 <            call calc_SC_prepair_rho(i, j, d, r, rijsq )
1368 >            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1369      endif
1370      
1371    end subroutine do_prepair
# Line 1338 | Line 1394 | contains
1394      real( kind = dp ) :: d(3), scaled(3)
1395      integer i
1396  
1397 <    d(1:3) = q_j(1:3) - q_i(1:3)
1397 >    d(1) = q_j(1) - q_i(1)
1398 >    d(2) = q_j(2) - q_i(2)
1399 >    d(3) = q_j(3) - q_i(3)
1400  
1401      ! Wrap back into periodic box if necessary
1402      if ( SIM_uses_PBC ) then
# Line 1361 | Line 1419 | contains
1419         else
1420            ! calc the scaled coordinates.
1421  
1364 –          do i = 1, 3
1365 –             scaled(i) = d(i) * HmatInv(i,i)
1422  
1423 <             ! wrap the scaled coordinates
1423 >          scaled(1) = d(1) * HmatInv(1,1)
1424 >          scaled(2) = d(2) * HmatInv(2,2)
1425 >          scaled(3) = d(3) * HmatInv(3,3)
1426 >          
1427 >          ! wrap the scaled coordinates
1428 >          
1429 >          scaled(1) = scaled(1) - dnint(scaled(1))
1430 >          scaled(2) = scaled(2) - dnint(scaled(2))
1431 >          scaled(3) = scaled(3) - dnint(scaled(3))
1432  
1433 <             scaled(i) = scaled(i) - anint(scaled(i))
1433 >          ! calc the wrapped real coordinates from the wrapped scaled
1434 >          ! coordinates
1435  
1436 <             ! calc the wrapped real coordinates from the wrapped scaled
1437 <             ! coordinates
1436 >          d(1) = scaled(1)*Hmat(1,1)
1437 >          d(2) = scaled(2)*Hmat(2,2)
1438 >          d(3) = scaled(3)*Hmat(3,3)
1439  
1374 –             d(i) = scaled(i)*Hmat(i,i)
1375 –          enddo
1440         endif
1441  
1442      endif
1443  
1444 <    r_sq = dot_product(d,d)
1444 >    r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1445  
1446    end subroutine get_interatomic_vector
1447  

Diff Legend

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