Skip to content

Commit

Permalink
Merge pull request #182 from openmopac/jimmy_pm6_org_update
Browse files Browse the repository at this point in the history 
Minor updates from Jimmy
  • Loading branch information
@godotalgorithm
godotalgorithm authored Sep 25, 2023
2 parents 2da75f6 + 67bad0e commit 23dff37
Show file tree
Hide file tree
Showing 8 changed files with 48 additions and 28 deletions.
4 changes: 2 additions & 2 deletions src/compfg.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -359,9 +359,9 @@ subroutine compfg(xparam, int, escf, fulscf, grad, lgrad)
hpress = 0.d0
if (Abs (pressure) > 1.d-4) then
if (id == 1) then
hpress = -pressure * Sqrt (dot(tvec(1, 1), tvec(1, 1), 3))
hpress = pressure * Sqrt (dot(tvec(1, 1), tvec(1, 1), 3))
else if (id == 3) then
hpress = -pressure * volume (tvec, 3)
hpress = pressure * volume (tvec, 3)
end if
atheat = atheat + hpress
end if
Expand Down
25 changes: 17 additions & 8 deletions src/forces/deriv.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,11 +254,20 @@ subroutine deriv(geo, gradnt)
! 1N = J/M = 10**(-3)/4.184 kcal/M = 4.184*10**(-3)*10**(-10) kcal/Angstrom
press = pressure / Sqrt (dot(tvec(1, 1), tvec(1, 1), 3))
do j = 1, 3
dxyz(j + i) = dxyz(j + i) + tvec(j, 1) * press
dxyz(j + i + 3) = dxyz(j + i + 3) - tvec(j, 1) * press
dxyz(j + i) = dxyz(j + i) - tvec(j, 1) * press
dxyz(j + i + 3) = dxyz(j + i + 3) + tvec(j, 1) * press
end do
else if (id == 3) then
! Transition vector derivatives of pressure times volume
!
! Derivatives are calculated in a redundant coordinate system where the
! atomic coordinates in the central cell plus translation vectors are
! replaced by atomic coordinates in the central cell plus atomic coordinates
! in all coupled non-central cells. The formula for volume that is being
! differentiated uses the coordinate difference between a central atom and
! its neighbors as proxies for the translation vectors in the volume formula.
! The choice of atom in this formula and the resulting derivatives are completely
! arbitrary and do not change the derivatives in the original coordinate system.
press = ((tvec(2, 1)*tvec(3, 2)-tvec(3, 1)*tvec(2, 2))*tvec(1, 3) + &
(tvec(3, 1)*tvec(1, 2)-tvec(1, 1)*tvec(3, 2))*tvec(2, 3) + &
(tvec(1, 1)*tvec(2, 2)-tvec(2, 1)*tvec(1, 2))*tvec(3, 3))
Expand All @@ -279,30 +288,30 @@ subroutine deriv(geo, gradnt)
!
i = 3*(l1u * (2*l2u+1) * (2*l3u+1) + l2u * (2*l3u+1) + l3u - 1)
do j = 1, 3
dxyz(j + i) = dxyz(j + i) + tderiv(j, 1)
dxyz(j + i) = dxyz(j + i) + tderiv(j, 2)
dxyz(j + i) = dxyz(j + i) + tderiv(j, 3)
dxyz(j + i) = dxyz(j + i) - tderiv(j, 1)
dxyz(j + i) = dxyz(j + i) - tderiv(j, 2)
dxyz(j + i) = dxyz(j + i) - tderiv(j, 3)
end do
!
! Cell in 0,0,1 position
!
i = 3*(l1u * (2*l2u+1) * (2*l3u+1) + l2u * (2*l3u+1) + l3u)
do j = 1, 3
dxyz(j + i) = dxyz(j + i) - tderiv(j, 3)
dxyz(j + i) = dxyz(j + i) + tderiv(j, 3)
end do
!
! Cell in 0,1,0 position
!
i = 3*(l1u * (2*l2u+1) * (2*l3u+1) + (l2u+1) * (2*l3u+1) + l3u - 1)
do j = 1, 3
dxyz(j + i) = dxyz(j + i) - tderiv(j, 2)
dxyz(j + i) = dxyz(j + i) + tderiv(j, 2)
end do
!
! Cell in 1,0,0 position
!
i = 3*((l1u+1) * (2*l2u+1) * (2*l3u+1) + l2u * (2*l3u+1) + l3u - 1)
do j = 1, 3
dxyz(j + i) = dxyz(j + i) - tderiv(j, 1)
dxyz(j + i) = dxyz(j + i) + tderiv(j, 1)
end do
end if
end if
Expand Down
7 changes: 4 additions & 3 deletions src/input/getdat.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ subroutine getdat(input, output)
!-----------------------------------------------
integer, parameter :: from_data_set = 7
integer :: i, j, io_stat, l, nlines
logical :: exists, arc_file, comments = .true.
logical :: exists, arc_file, comments = .true., double_plus
character :: line1*3000, num1*1, num2*1
character, allocatable :: tmp_comments(:)*120
double precision, external :: reada
Expand Down Expand Up @@ -294,7 +294,8 @@ subroutine getdat(input, output)
i = index(keywrd, "GEO-DAT")
if (i /= 0) keywrd(i:i+6) = "GEO_DAT"
i = index(keywrd, "GEO-REF")
if (i /= 0) keywrd(i:i+6) = "GEO_REF"
if (i /= 0) keywrd(i:i+6) = "GEO_REF"
double_plus = (index(keywrd, " ++ ") /= 0)
if (index(keywrd, " GEO_DAT") + index(keywrd, " SETUP")/= 0) then
nlines = nlines + 3
else if (.not. is_PARAM .and. nlines < 4) then
Expand Down Expand Up @@ -382,7 +383,7 @@ subroutine getdat(input, output)
line = ' '
write (input, '(A241)') line
rewind input
1000 if (nlines < 3 .and. .not. is_PARAM) then
1000 if (nlines < 3 .and. .not. is_PARAM .and. line1 == " " .and. .not. double_plus) then
inquire(unit=output, opened=exists)
if (.not. exists) open(unit=output, file=trim(jobnam)//'.out')
#ifdef MOPAC_F2003
Expand Down
15 changes: 13 additions & 2 deletions src/input/wrtkey.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1168,6 +1168,17 @@ subroutine wrtcon (allkey)
end do
end if
line = trim(allkey)
i = 0
if (myword(line, " OPT-")) i = i + 1
if (myword(line, " OPT ")) i = i + 1
if (myword(line, " OPT(")) i = i + 1
if (myword(line, " OPT=")) i = i + 1
if (i > 1) then
call mopend("MORE THAN ONE TYPE OF ""OPT"" REQUESTED. THIS IS NOT ALLOWED")
write(iw,'(/10x,a,//,a,/)')"KEYWORDS SUPPLIED:", trim(keywrd)
return
end if
line = trim(allkey)
if (myword(allkey, " OPT-") .or. myword(allkey, " OPT ") .or. myword(allkey, " OPT(") &
.or. myword(allkey, " OPT=")) then
i = index(line, " OPT=")
Expand Down Expand Up @@ -1240,7 +1251,7 @@ subroutine wrtcon (allkey)
pressure = reada (keywrd, Index (keywrd, " P="))
if (id == 1) then
write (iw, '(" * P= - TENSION IN POLYMER=", g13.6, " NEWTONS PER MOLE")') pressure
pressure = -pressure * 10.d0 ** (-13) / 4.184d0
pressure = pressure * 10.d0 ** (-13) / 4.184d0
else if (id == 2) then
else if (id == 3) then
i = Index (keywrd, " P=")
Expand All @@ -1259,7 +1270,7 @@ subroutine wrtcon (allkey)
! Divide by 4184 to convert from J/M**3/mol to Kcal/M**3/mol
! Divide by 10**30 to convert from Kcal/M**3/mol to Kcal/Angstrom**3/mol
!
pressure = -(fpcref(1,10)*pressure) / (4184.d0*10.d0**30)
pressure = (fpcref(1,10)*pressure) / (4184.d0*10.d0**30)
else
write (iw, *) " Keyword 'P=n.nn' is not allowed here"
call mopend("Keyword 'P=n.nn' is not allowed here")
Expand Down
5 changes: 2 additions & 3 deletions src/models/refer.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,9 +148,8 @@ subroutine refer
end if
write (iw, '(/,a)') " General Reference for PM6-ORG:"
write (iw, '(1x,a)') &
& """A semiempirical method optimized for modeling proteins""", &
& "J. J. P. Stewart and A. C. Stewart, J. Mol. Model. 29, 284 (2023) ", &
& "https://doi.org/10.1007/s00894-023-05695-1"
& """A semiempirical method optimized for modeling proteins"", J. J. P. Stewart and A. C. Stewart", &
& "Journal of Molecular Modeling (2023) 29:284 [https://doi.org/10.1007/s00894-023-05695-1]"
return
end if
mixok = index(keywrd,'PARASOK') /= 0
Expand Down
6 changes: 3 additions & 3 deletions src/moldat.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1308,9 +1308,9 @@ subroutine calculate_voigt
voigt(6) = voigt(6) + 0.5*(dsum(2)*dsum1(1) + dsum(1)*dsum1(2))
end if
end do
! Convert from enthalpy to internal energy (sign is flipped on pressure)
! Convert from enthalpy to internal energy
do i = 1, 3
voigt(i) = voigt(i) + pressure*volume (tvec, 3)
voigt(i) = voigt(i) - pressure*volume (tvec, 3)
end do
! Convert units of stress tensor to GPa
xi = 1.d-9*(4184.d0*10.d0**30)/(fpc_10*volume (tvec, 3))
Expand Down Expand Up @@ -1379,7 +1379,7 @@ subroutine write_pressure(iprt)
end if
i2 = 1
end if
xi = (xi - pressure * (4184.d0*10.d0**30)/ fpc_10)*1.d-9
xi = (xi + pressure * (4184.d0*10.d0**30)/ fpc_10)*1.d-9
ndim = ndim + 1
press(ndim) = xi
write(line,'(10x,a,i4,a,f7.2,a)')"Tv(", k,") Pressure:",xi," GPa"
Expand Down
2 changes: 1 addition & 1 deletion src/molkst_C.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,7 @@ module molkst_C
& cupper, & ! Upper bound for transition to point charges in solid-state
! Default cutofp
!
& pressure, & ! Term Pressure or stress for solid-state and polymer work (sign flipped)
& pressure, & ! Term Pressure or stress for solid-state and polymer work
! Definition Applied isotropic pressure (for solids) or pull (for polymers)
! Units Pascals (Newtons per square meter) (solids) or Newtons per Mole (polymers)
! Default 0.0
Expand Down
12 changes: 6 additions & 6 deletions src/output/writmo.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,9 +162,9 @@ subroutine writmo
! Remove energy term arising from the external pressure,
!
if( id == 1) then
escf = escf + pressure * dsqrt(dot(tvec(1,1), tvec(1,1),3))
escf = escf - pressure * dsqrt(dot(tvec(1,1), tvec(1,1),3))
else if (id == 3) then
escf = escf + pressure * volume(tvec,3)
escf = escf - pressure * volume(tvec,3)
end if
end if
if (use_ref_geo) then
Expand Down Expand Up @@ -281,13 +281,13 @@ subroutine writmo
hpress = 0.d0
if (Abs (pressure) > 1.d-4) then
if (id == 1) then
hpress = -pressure * dsqrt (dot(tvec(1, 1), tvec(1, 1), 3))
hpress = pressure * dsqrt (dot(tvec(1, 1), tvec(1, 1), 3))
else if (id == 3) then
sum = volume (tvec, 3)
hpress = -pressure * sum
hpress = pressure * sum
write (iw, '( 10X,''ENERGY DUE TO PRESSURE ='',F17.5,'' KCAL/MOL'' )') &
hpress
sum = -(4184.d0*10.d0**30)*pressure/fpc_10
sum = (4184.d0*10.d0**30)*pressure/fpc_10
if (abs(sum) > 1.d9) then
write (iw, '(10X,''PRESSURE ='',F17.5,'' Gp'' )') sum*1.d-9
else
Expand Down Expand Up @@ -1081,7 +1081,7 @@ subroutine writmo
sum = volume (tvec, 3)
write (iwrite, '( 10X,''ENERGY DUE TO PRESSURE ='',F17.5,'' KCAL/MOL'' )') &
hpress
sum = -(4184.d0*10.d0**30)*pressure/fpc_10
sum = (4184.d0*10.d0**30)*pressure/fpc_10
if (abs(sum) > 1.d9) then
write (iwrite, '(10X,''PRESSURE ='',F17.5,'' Gp'' )') sum*1.d-9
else
Expand Down

0 comments on commit 23dff37

Please sign in to comment.