Skip to content

Commit

Permalink
Revert non-standard internal sign change of pressure
Browse files Browse the repository at this point in the history
  • Loading branch information
@godotalgorithm
godotalgorithm committed Sep 25, 2023
1 parent 65ad524 commit 67bad0e
Show file tree
Hide file tree
Showing 6 changed files with 22 additions and 22 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))

Check warning on line 362 in src/compfg.F90

View check run for this annotation

Codecov / codecov/patch

src/compfg.F90#L362 

Added line #L362 was not covered by tests
else if (id == 3) then
hpress = -pressure * volume (tvec, 3)
hpress = pressure * volume (tvec, 3)

Check warning on line 364 in src/compfg.F90

View check run for this annotation

Codecov / codecov/patch

src/compfg.F90#L364 

Added line #L364 was not covered by tests
end if
atheat = atheat + hpress
end if
Expand Down
16 changes: 8 additions & 8 deletions src/forces/deriv.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,8 +254,8 @@ 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

Check warning on line 258 in src/forces/deriv.F90

View check run for this annotation

Codecov / codecov/patch

src/forces/deriv.F90#L257-L258 

Added lines #L257 - L258 were not covered by tests
end do
else if (id == 3) then
! Transition vector derivatives of pressure times volume
Expand Down Expand Up @@ -288,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)

Check warning on line 293 in src/forces/deriv.F90

View check run for this annotation

Codecov / codecov/patch

src/forces/deriv.F90#L291-L293 

Added lines #L291 - L293 were not covered by tests
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)

Check warning on line 300 in src/forces/deriv.F90

View check run for this annotation

Codecov / codecov/patch

src/forces/deriv.F90#L300 

Added line #L300 was not covered by tests
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)

Check warning on line 307 in src/forces/deriv.F90

View check run for this annotation

Codecov / codecov/patch

src/forces/deriv.F90#L307 

Added line #L307 was not covered by tests
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)

Check warning on line 314 in src/forces/deriv.F90

View check run for this annotation

Codecov / codecov/patch

src/forces/deriv.F90#L314 

Added line #L314 was not covered by tests
end do
end if
end if
Expand Down
4 changes: 2 additions & 2 deletions src/input/wrtkey.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1251,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

Check warning on line 1254 in src/input/wrtkey.F90

View check run for this annotation

Codecov / codecov/patch

src/input/wrtkey.F90#L1254 

Added line #L1254 was not covered by tests
else if (id == 2) then
else if (id == 3) then
i = Index (keywrd, " P=")
Expand All @@ -1270,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)

Check warning on line 1273 in src/input/wrtkey.F90

View check run for this annotation

Codecov / codecov/patch

src/input/wrtkey.F90#L1273 

Added line #L1273 was not covered by tests
else
write (iw, *) " Keyword 'P=n.nn' is not allowed here"
call mopend("Keyword 'P=n.nn' is not allowed here")
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)

Check warning on line 1313 in src/moldat.F90

View check run for this annotation

Codecov / codecov/patch

src/moldat.F90#L1313 

Added line #L1313 was not covered by tests
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

Check warning on line 1382 in src/moldat.F90

View check run for this annotation

Codecov / codecov/patch

src/moldat.F90#L1382 

Added line #L1382 was not covered by tests
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))

Check warning on line 165 in src/output/writmo.F90

View check run for this annotation

Codecov / codecov/patch

src/output/writmo.F90#L165 

Added line #L165 was not covered by tests
else if (id == 3) then
escf = escf + pressure * volume(tvec,3)
escf = escf - pressure * volume(tvec,3)

Check warning on line 167 in src/output/writmo.F90

View check run for this annotation

Codecov / codecov/patch

src/output/writmo.F90#L167 

Added line #L167 was not covered by tests
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))

Check warning on line 284 in src/output/writmo.F90

View check run for this annotation

Codecov / codecov/patch

src/output/writmo.F90#L284 

Added line #L284 was not covered by tests
else if (id == 3) then
sum = volume (tvec, 3)
hpress = -pressure * sum
hpress = pressure * sum

Check warning on line 287 in src/output/writmo.F90

View check run for this annotation

Codecov / codecov/patch

src/output/writmo.F90#L287 

Added line #L287 was not covered by tests
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

Check warning on line 290 in src/output/writmo.F90

View check run for this annotation

Codecov / codecov/patch

src/output/writmo.F90#L290 

Added line #L290 was not covered by tests
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

Check warning on line 1084 in src/output/writmo.F90

View check run for this annotation

Codecov / codecov/patch

src/output/writmo.F90#L1084 

Added line #L1084 was not covered by tests
if (abs(sum) > 1.d9) then
write (iwrite, '(10X,''PRESSURE ='',F17.5,'' Gp'' )') sum*1.d-9
else
Expand Down

0 comments on commit 67bad0e

Please sign in to comment.