45 |
|
|
46 |
|
!! @author Charles F. Vardeman II |
47 |
|
!! @author Matthew Meineke |
48 |
< |
!! @version $Id: doForces.F90,v 1.71 2005-12-15 21:43:16 gezelter Exp $, $Date: 2005-12-15 21:43:16 $, $Name: not supported by cvs2svn $, $Revision: 1.71 $ |
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 |
64 |
|
use eam |
65 |
|
use suttonchen |
66 |
|
use status |
67 |
+ |
use interpolation |
68 |
|
#ifdef IS_MPI |
69 |
|
use mpiSimulation |
70 |
|
#endif |
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 |
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 |
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 |
|
|
354 |
|
if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then |
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)) |
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 |
585 |
|
localError = 0 |
586 |
|
call setLJDefaultCutoff( defaultRcut, defaultDoShift ) |
587 |
|
call setElectrostaticCutoffRadius( defaultRcut, defaultRsw ) |
588 |
< |
call setCutoffEAM( defaultRcut, localError) |
589 |
< |
if (localError /= 0) then |
584 |
< |
write(errMsg, *) 'An error has occured in setting the EAM cutoff' |
585 |
< |
call handleError("setCutoffs", errMsg) |
586 |
< |
end if |
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 |
|
haveGtypeCutoffMap = .false. |
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 |
|
|
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 |
|
|
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) |
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. |
746 |
|
|
747 |
|
if (FF_uses_EAM) then |
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, & |
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 |
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 |
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 |
1419 |
|
else |
1420 |
|
! calc the scaled coordinates. |
1421 |
|
|
1400 |
– |
do i = 1, 3 |
1401 |
– |
scaled(i) = d(i) * HmatInv(i,i) |
1422 |
|
|
1423 |
< |
! wrap the scaled coordinates |
1424 |
< |
|
1425 |
< |
scaled(i) = scaled(i) - anint(scaled(i)) |
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 |
< |
! calc the wrapped real coordinates from the wrapped scaled |
1434 |
< |
! coordinates |
1433 |
> |
! calc the wrapped real coordinates from the wrapped scaled |
1434 |
> |
! coordinates |
1435 |
|
|
1436 |
< |
d(i) = scaled(i)*Hmat(i,i) |
1437 |
< |
enddo |
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 |
> |
|
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 |
|
|