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 2592 by gezelter, Thu Feb 16 21:40:20 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.76 2006-02-16 21:40:20 gezelter Exp $, $Date: 2006-02-16 21:40:20 $, $Name: not supported by cvs2svn $, $Revision: 1.76 $
49  
50  
51   module doForces
# Line 87 | Line 87 | module doForces
87    logical, save :: haveInteractionHash = .false.
88    logical, save :: haveGtypeCutoffMap = .false.
89    logical, save :: haveDefaultCutoffs = .false.
90 <  logical, save :: haveRlist = .false.
90 >  logical, save :: haveSkinThickness = .false.
91 >  logical, save :: haveElectrostaticSummationMethod = .false.
92 >  logical, save :: haveCutoffPolicy = .false.
93 >  logical, save :: VisitCutoffsAfterComputing = .false.
94  
95    logical, save :: FF_uses_DirectionalAtoms
96    logical, save :: FF_uses_Dipoles
# Line 106 | Line 109 | module doForces
109    logical, save :: SIM_uses_PBC
110  
111    integer, save :: electrostaticSummationMethod
112 +  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
113  
114 +  real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
115 +  real(kind=dp), save :: skinThickness
116 +  logical, save :: defaultDoShift
117 +
118    public :: init_FF
119 <  public :: setDefaultCutoffs
119 >  public :: setCutoffs
120 >  public :: cWasLame
121 >  public :: setElectrostaticMethod
122 >  public :: setCutoffPolicy
123 >  public :: setSkinThickness
124    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
125  
126   #ifdef PROFILE
127    public :: getforcetime
# Line 144 | Line 149 | module doForces
149    end type gtypeCutoffs
150    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
151  
147 –  integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
148 –  real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
149 –  real(kind=dp),save :: listSkin
150 –  
152   contains
153  
154 <  subroutine createInteractionHash(status)
154 >  subroutine createInteractionHash()
155      integer :: nAtypes
155 –    integer, intent(out) :: status
156      integer :: i
157      integer :: j
158      integer :: iHash
# Line 176 | Line 176 | contains
176      logical :: j_is_SC
177      logical :: j_is_MEAM
178      real(kind=dp) :: myRcut
179 –
180 –
181 –    status = 0  
179  
180      if (.not. associated(atypes)) then
181 <       call handleError("atype", "atypes was not present before call of createInteractionHash!")
185 <       status = -1
181 >       call handleError("doForces", "atypes was not present before call of createInteractionHash!")
182         return
183      endif
184      
185      nAtypes = getSize(atypes)
186      
187      if (nAtypes == 0) then
188 <       status = -1
188 >       call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
189         return
190      end if
191  
# Line 276 | Line 272 | contains
272      haveInteractionHash = .true.
273    end subroutine createInteractionHash
274  
275 <  subroutine createGtypeCutoffMap(stat)
275 >  subroutine createGtypeCutoffMap()
276  
281 –    integer, intent(out), optional :: stat
277      logical :: i_is_LJ
278      logical :: i_is_Elect
279      logical :: i_is_Sticky
# Line 286 | Line 281 | contains
281      logical :: i_is_GB
282      logical :: i_is_EAM
283      logical :: i_is_Shape
284 +    logical :: i_is_SC
285      logical :: GtypeFound
286  
287      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
# Line 293 | Line 289 | contains
289      integer :: nGroupsInRow
290      integer :: nGroupsInCol
291      integer :: nGroupTypesRow,nGroupTypesCol
292 <    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
292 >    real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
293      real(kind=dp) :: biggestAtypeCutoff
294  
299 –    stat = 0
295      if (.not. haveInteractionHash) then
296 <       call createInteractionHash(myStatus)      
302 <       if (myStatus .ne. 0) then
303 <          write(default_error, *) 'createInteractionHash failed in doForces!'
304 <          stat = -1
305 <          return
306 <       endif
296 >       call createInteractionHash()      
297      endif
298   #ifdef IS_MPI
299      nGroupsInRow = getNgroupsInRow(plan_group_row)
# Line 321 | Line 311 | contains
311            call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
312            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
313            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
314 <          
314 >          call getElementProperty(atypes, i, "is_SC", i_is_SC)
315  
316            if (haveDefaultCutoffs) then
317               atypeMaxCutoff(i) = defaultRcut
# Line 354 | Line 344 | contains
344                  thisRcut = getShapeCut(i)
345                  if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
346               endif
347 +             if (i_is_SC) then
348 +                thisRcut = getSCCut(i)
349 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
350 +             endif
351            endif
352 <          
359 <          
352 >                    
353            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
354               biggestAtypeCutoff = atypeMaxCutoff(i)
355            endif
356  
357         endif
358      enddo
366 –  
367 –
359      
360      istart = 1
361      jstart = 1
# Line 408 | Line 399 | contains
399         allocate(groupToGtypeCol(jend))
400      end if
401  
402 <    if(.not.associated(groupToGtypeCol)) then
403 <       allocate(groupToGtypeCol(jend))
402 >    if(.not.associated(groupMaxCutoffCol)) then
403 >       allocate(groupMaxCutoffCol(jend))
404      else
405 <       deallocate(groupToGtypeCol)
406 <       allocate(groupToGtypeCol(jend))
405 >       deallocate(groupMaxCutoffCol)
406 >       allocate(groupMaxCutoffCol(jend))
407      end if
408      if(.not.associated(gtypeMaxCutoffCol)) then
409         allocate(gtypeMaxCutoffCol(jend))
# Line 435 | Line 426 | contains
426      
427      tol = 1.0d-6
428      nGroupTypesRow = 0
429 <
429 >    nGroupTypesCol = 0
430      do i = istart, iend      
431         n_in_i = groupStartRow(i+1) - groupStartRow(i)
432         groupMaxCutoffRow(i) = 0.0_dp
# Line 450 | Line 441 | contains
441               groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
442            endif          
443         enddo
453 –
444         if (nGroupTypesRow.eq.0) then
445            nGroupTypesRow = nGroupTypesRow + 1
446            gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
# Line 513 | Line 503 | contains
503      groupMaxCutoffCol => groupMaxCutoffRow
504   #endif
505  
516 –
517 –
518 –
519 –
506      !! allocate the gtypeCutoffMap here.
507      allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
508      !! then we do a double loop over all the group TYPES to find the cutoff
509      !! map between groups of two types
510      tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
511  
512 <    do i = 1, nGroupTypesRow
512 >    do i = 1, nGroupTypesRow      
513         do j = 1, nGroupTypesCol
514        
515            select case(cutoffPolicy)
# Line 538 | Line 524 | contains
524               return
525            end select
526            gtypeCutoffMap(i,j)%rcut = thisRcut
527 +          
528 +          if (thisRcut.gt.largestRcut) largestRcut = thisRcut
529 +
530            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
531  
532 +          if (.not.haveSkinThickness) then
533 +             skinThickness = 1.0_dp
534 +          endif
535 +
536 +          gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
537 +
538            ! sanity check
539  
540            if (haveDefaultCutoffs) then
# Line 552 | Line 544 | contains
544            endif
545         enddo
546      enddo
547 +
548      if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
549      if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
550      if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
# Line 565 | Line 558 | contains
558      haveGtypeCutoffMap = .true.
559     end subroutine createGtypeCutoffMap
560  
561 <   subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
569 <     real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
570 <     integer, intent(in) :: cutPolicy
561 >   subroutine setCutoffs(defRcut, defRsw)
562  
563 +     real(kind=dp),intent(in) :: defRcut, defRsw
564 +     character(len = statusMsgSize) :: errMsg
565 +     integer :: localError
566 +
567       defaultRcut = defRcut
568       defaultRsw = defRsw
569 <     defaultRlist = defRlist
570 <     cutoffPolicy = cutPolicy
569 >    
570 >     defaultDoShift = .false.
571 >     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
572 >        
573 >        write(errMsg, *) &
574 >             'cutoffRadius and switchingRadius are set to the same', newline &
575 >             // tab, 'value.  OOPSE will use shifted ', newline &
576 >             // tab, 'potentials instead of switching functions.'
577 >        
578 >        call handleInfo("setCutoffs", errMsg)
579 >        
580 >        defaultDoShift = .true.
581 >        
582 >     endif
583  
584 +     localError = 0
585 +     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
586 +     call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
587 +     call setCutoffEAM( defaultRcut, localError)
588 +     if (localError /= 0) then
589 +       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
590 +       call handleError("setCutoffs", errMsg)
591 +     end if
592 +     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
593 +     call setHmatDangerousRcutValue(defaultRcut)
594 +
595       haveDefaultCutoffs = .true.
596 <   end subroutine setDefaultCutoffs
596 >     haveGtypeCutoffMap = .false.
597 >   end subroutine setCutoffs
598  
599 +   subroutine cWasLame()
600 +    
601 +     VisitCutoffsAfterComputing = .true.
602 +     return
603 +    
604 +   end subroutine cWasLame
605 +  
606     subroutine setCutoffPolicy(cutPolicy)
607 <
607 >    
608       integer, intent(in) :: cutPolicy
609 +    
610       cutoffPolicy = cutPolicy
611 <     call createGtypeCutoffMap()
612 <   end subroutine setCutoffPolicy
586 <    
611 >     haveCutoffPolicy = .true.
612 >     haveGtypeCutoffMap = .false.
613      
614 <  subroutine setSimVariables()
615 <    SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
616 <    SIM_uses_EAM = SimUsesEAM()
591 <    SIM_uses_SC  = SimUsesSC()
592 <    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
593 <    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
594 <    SIM_uses_PBC = SimUsesPBC()
614 >   end subroutine setCutoffPolicy
615 >  
616 >   subroutine setElectrostaticMethod( thisESM )
617  
618 <    haveSIMvariables = .true.
618 >     integer, intent(in) :: thisESM
619  
620 <    return
621 <  end subroutine setSimVariables
620 >     electrostaticSummationMethod = thisESM
621 >     haveElectrostaticSummationMethod = .true.
622 >    
623 >   end subroutine setElectrostaticMethod
624  
625 +   subroutine setSkinThickness( thisSkin )
626 +    
627 +     real(kind=dp), intent(in) :: thisSkin
628 +    
629 +     skinThickness = thisSkin
630 +     haveSkinThickness = .true.    
631 +     haveGtypeCutoffMap = .false.
632 +    
633 +   end subroutine setSkinThickness
634 +      
635 +   subroutine setSimVariables()
636 +     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
637 +     SIM_uses_EAM = SimUsesEAM()
638 +     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
639 +     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
640 +     SIM_uses_PBC = SimUsesPBC()
641 +     SIM_uses_SC = SimUsesSC()
642 +    
643 +     haveSIMvariables = .true.
644 +    
645 +     return
646 +   end subroutine setSimVariables
647 +
648    subroutine doReadyCheck(error)
649      integer, intent(out) :: error
650  
# Line 606 | Line 653 | contains
653      error = 0
654  
655      if (.not. haveInteractionHash) then      
656 <       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
656 >       call createInteractionHash()      
657      endif
658  
659      if (.not. haveGtypeCutoffMap) then        
660 <       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
660 >       call createGtypeCutoffMap()      
661      endif
662  
663 +
664 +    if (VisitCutoffsAfterComputing) then
665 +       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
666 +       call setHmatDangerousRcutValue(largestRcut)
667 +    endif
668 +
669 +
670      if (.not. haveSIMvariables) then
671         call setSimVariables()
672      endif
# Line 658 | Line 700 | contains
700    end subroutine doReadyCheck
701  
702  
703 <  subroutine init_FF(thisESM, thisStat)
703 >  subroutine init_FF(thisStat)
704  
663 –    integer, intent(in) :: thisESM
705      integer, intent(out) :: thisStat  
706      integer :: my_status, nMatches
707      integer, pointer :: MatchList(:) => null()
# Line 668 | Line 709 | contains
709      !! assume things are copacetic, unless they aren't
710      thisStat = 0
711  
671 –    electrostaticSummationMethod = thisESM
672 –
712      !! init_FF is called *after* all of the atom types have been
713      !! defined in atype_module using the new_atype subroutine.
714      !!
# Line 680 | Line 719 | contains
719      FF_uses_Dipoles = .false.
720      FF_uses_GayBerne = .false.
721      FF_uses_EAM = .false.
722 +    FF_uses_SC = .false.
723  
724      call getMatchingElementList(atypes, "is_Directional", .true., &
725           nMatches, MatchList)
# Line 696 | Line 736 | contains
736      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
737      if (nMatches .gt. 0) FF_uses_EAM = .true.
738  
739 +    call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
740 +    if (nMatches .gt. 0) FF_uses_SC = .true.
741  
742 +
743      haveSaneForceField = .true.
744  
745      if (FF_uses_EAM) then
# Line 766 | Line 809 | contains
809      real( kind = DP ) :: rVal
810      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
811      real(kind=dp) :: rfpot, mu_i, virial
812 +    real(kind=dp):: rCut
813      integer :: me_i, me_j, n_in_i, n_in_j
814      logical :: is_dp_i
815      integer :: neighborListSize
# Line 839 | Line 883 | contains
883         ! (but only on the first time through):
884         if (loop .eq. loopStart) then
885   #ifdef IS_MPI
886 <          call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
886 >          call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
887                 update_nlist)
888   #else
889 <          call checkNeighborList(nGroups, q_group, listSkin, &
889 >          call checkNeighborList(nGroups, q_group, skinThickness, &
890                 update_nlist)
891   #endif
892         endif
# Line 922 | Line 966 | contains
966  
967                     list(nlist) = j
968                  endif
969 +
970 +
971                  
972                  if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
973  
974 +                   rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
975                     if (loop .eq. PAIR_LOOP) then
976                        vij = 0.0d0
977                        fij(1:3) = 0.0d0
978                     endif
979                    
980 <                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
981 <                        in_switching_region)
980 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
981 >                        group_switch, in_switching_region)
982                    
983                     n_in_j = groupStartCol(j+1) - groupStartCol(j)
984                    
# Line 961 | Line 1008 | contains
1008                           if (loop .eq. PREPAIR_LOOP) then
1009   #ifdef IS_MPI                      
1010                              call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1011 <                                 rgrpsq, d_grp, do_pot, do_stress, &
1011 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1012                                   eFrame, A, f, t, pot_local)
1013   #else
1014                              call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1015 <                                 rgrpsq, d_grp, do_pot, do_stress, &
1015 >                                 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1016                                   eFrame, A, f, t, pot)
1017   #endif                                              
1018                           else
1019   #ifdef IS_MPI                      
1020                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1021                                   do_pot, eFrame, A, f, t, pot_local, vpair, &
1022 <                                 fpair, d_grp, rgrp)
1022 >                                 fpair, d_grp, rgrp, rCut)
1023   #else
1024                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1025                                   do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1026 <                                 d_grp, rgrp)
1026 >                                 d_grp, rgrp, rCut)
1027   #endif
1028                              vij = vij + vpair
1029                              fij(1:3) = fij(1:3) + fpair(1:3)
# Line 1184 | Line 1231 | contains
1231    end subroutine do_force_loop
1232  
1233    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1234 <       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp)
1234 >       eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1235  
1236      real( kind = dp ) :: vpair, sw
1237      real( kind = dp ), dimension(LR_POT_TYPES) :: pot
# Line 1201 | Line 1248 | contains
1248      real ( kind = dp ), intent(inout) :: r_grp
1249      real ( kind = dp ), intent(inout) :: d(3)
1250      real ( kind = dp ), intent(inout) :: d_grp(3)
1251 +    real ( kind = dp ), intent(inout) :: rCut
1252      real ( kind = dp ) :: r
1253      integer :: me_i, me_j
1254  
# Line 1221 | Line 1269 | contains
1269      iHash = InteractionHash(me_i, me_j)
1270      
1271      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1272 <       call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1272 >       call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1273              pot(VDW_POT), f, do_pot)
1274      endif
1275      
1276      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1277 <       call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1277 >       call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1278              pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1279      endif
1280      
# Line 1246 | Line 1294 | contains
1294      endif
1295      
1296      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1297 <       call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1297 >       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1298              pot(VDW_POT), A, f, t, do_pot)
1299      endif
1300      
# Line 1266 | Line 1314 | contains
1314      endif
1315  
1316      if ( iand(iHash, SC_PAIR).ne.0 ) then      
1317 <       call do_SC_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1317 >       call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1318              pot(METALLIC_POT), f, do_pot)
1319      endif
1320  
# Line 1274 | Line 1322 | contains
1322      
1323    end subroutine do_pair
1324  
1325 <  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1325 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1326         do_pot, do_stress, eFrame, A, f, t, pot)
1327  
1328      real( kind = dp ) :: sw
# Line 1286 | Line 1334 | contains
1334  
1335      logical, intent(inout) :: do_pot, do_stress
1336      integer, intent(in) :: i, j
1337 <    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1337 >    real ( kind = dp ), intent(inout)    :: rijsq, rcijsq, rCut
1338      real ( kind = dp )                :: r, rc
1339      real ( kind = dp ), intent(inout) :: d(3), dc(3)
1340  
# Line 1305 | Line 1353 | contains
1353      iHash = InteractionHash(me_i, me_j)
1354  
1355      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1356 <            call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1356 >            call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1357      endif
1358  
1359      if ( iand(iHash, SC_PAIR).ne.0 ) then      
1360 <            call calc_SC_prepair_rho(i, j, d, r, rijsq )
1360 >            call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1361      endif
1362      
1363    end subroutine do_prepair

Diff Legend

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