Page 1 of 2 1 2 >
Topic Options
#25199 - 08/27/10 03:23 AM Pressure unit in charmm
Giselle Offline
Forum Member

Registered: 08/03/10
Posts: 8
Hi everyone,

I've got a question about the unit of pressure in charmm. I did a simulation of a box of water molecules in NPT ensemble by using 'cpt' command lines as below:

pcons pmass 360 pgamma 20 pref 1.0 - ! since there are 1000 water molecules
! in the system, I set pmass as 360

I am sure the 'pref' value uses the unit of 'atm'. But when I read my output file to check if the simulation result was reasonable, I was just confused by the unit it used.

In the output file, the system pressure is recorded like this:

AVER PRESS: VIRE VIRI PRESSE PRESSI VOLUme
LAVE PRESS> 478.01444 -1155.66231 -1112.90548 23.11227 29451.39517

The value of inert pressure 23.11227 is definitely not in the unit of atm. I tried to convert it into an AKMA unit, but it is still not reasonable.

1 atm = 101325 N/m^2
1 N = 1.4388 kcal/mol-A
1 m^2 = 1E20 A^2

1 atm = 101325*1.4388*1E-20 kcal/mol-A^3 =1.4579E-15 kcal/mol-A^3

Are there any mistakes in my unit conversion? Or just because in the output file it is a totally different unit for pressure?
If it really uses atm for pressure in output file, then what was wrong with my simulation?

I thought it would be a very basic question, but I just couldn't find any relevant topic after searching through the forum.

Thanks for any suggestions!

Giselle

Top
#25200 - 08/27/10 03:52 AM Re: Pressure unit in charmm [Re: Giselle]
lennart Offline

Forum Member

Registered: 09/25/03
Posts: 4688
Loc: ~ 59N, 15E
It is atm. Why do you think it isn't? If you expected 1.000 it may be that your simulation is too short. Note that for all practical intents and purposes 1 atm is identical to 23 atm, which is the pressure at 230 m depth in the sea, and I do not think organisms living in this depth range have any differences in their proteins due to the pressure. Pressure experiments on proteins need thousands of atmospheres to get any effect at all.
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#25201 - 08/27/10 05:35 AM Re: Pressure unit in charmm [Re: lennart]
Giselle Offline
Forum Member

Registered: 08/03/10
Posts: 8
Hi Lennart,

Thanks for your replying!

And about my simulation, I separated it into three stages. First was a 1 ns long simulation in NVT-hoover ensemble, and then restarted it as a 500 ps simulation in NPT-hoover ensemble, the above two stages were used to equilibrate the system. After that, I ran a 5 ns simulation to get the trajectory by restarting the previous simulation, which was the third stage.

Actually, I got this average internal pressure (23.11227 atm) after the whole 6.5 ns simulation. The average internal pressure was -79 atm right after the second stage. So, is it because the fluctuation of internal pressure being so big ( about 400 ) that there's no need to expect the value to be exactly 1 atm ? Or simply because 6.5 ns is not long enough to reach 1 atm?

Thanks for further suggestions.

Giselle

Top
#25203 - 08/27/10 05:49 AM Re: Pressure unit in charmm [Re: Giselle]
lennart Offline

Forum Member

Registered: 09/25/03
Posts: 4688
Loc: ~ 59N, 15E
6.5 ns is probably long enough. Did you actually calculate the average over 6.5ns? CHARMM prints the number of steps over which the average was calculated just before the LAVE> lines.
You may have to calculate the real long-term average yourself from the data in the logfile; there are examples in the Script Archive that show how to do this.
_________________________
Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden

Top
#25204 - 08/27/10 06:21 AM Re: Pressure unit in charmm [Re: lennart]
Giselle Offline
Forum Member

Registered: 08/03/10
Posts: 8
Hi Lennart,

Oh, I did just look at the LAVE> part, which only averaged the last 5000 steps of the simulation since I thought it was good enough to check my simulation result...

Well, after I tried averaging the pressure from the data in the output file, it did become smaller, but is still not exactly 1 atm.

Stage__Time (ns)____Ensemble used____Average Pressure (atm)
1______0-1_________NVT_______________17027_______________
2______1-1.5_______NPT__________________12.75_____________
3______1.5-6.5_____NVT__________________62.53_____________

However, I think that's because the way I arranged these three stages needs some corrections. The first stage had kind of exploded...

Thanks for your advices! I'll do another try!

Giselle



Edited by Giselle (08/27/10 06:22 AM)

Top
#25211 - 08/27/10 05:53 PM Re: Pressure unit in charmm [Re: Giselle]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8259
Loc: 39 03 48 N, 77 06 54 W
It's usually a good idea to start with NPT, and run a couple ns to allow pressure equilibration, before imposing a constant volume.

I recently posted a suggested method for changing from NPT to NVT here
_________________________
Rick Venable
computational chemist


Top
#25217 - 08/28/10 01:12 AM Re: Pressure unit in charmm [Re: rmv]
Kenno Offline
Forum Member

Registered: 12/19/05
Posts: 1535
Loc: Baltimore, MD
Or one could just run the whole simulation with NPT. We always recommend beginners to do so - less possibility for things to go wrong. And the computational overhead of the barostat is not that dramatic...

Top
#37147 - 09/07/18 06:22 AM Re: Pressure unit in charmm [Re: Giselle]
vidhya Offline
Forum Member

Registered: 08/09/18
Posts: 10
Loc: India
Hi, I dont know if I should post my query in this but since this is related to pressure in my NPT simulation, I would think it is fine.

I am trying to run simulation using SSBP (Spherical Solvent Boundary Potential) and using HOOVER thermostat to maintain my temperature.
When I try to use PCONS command, it says pcons can be used only for CRYSTAL facility.
Since I am not using PBC, and SSBP by itself says that it is a constant temperature and pressure system, could you suggest how to define NPT for my system?

Currently my input file is like this.

!Defining spherical solvent boundary potential
MMFP
SSBP SELE TYPE OH2 END
GEO sphere RCM -
xref 0.0 yref 0.0 zref 0.0 -
force 10.0 droff 0.0 -
select segid PI .or. segid NI end
END

!Defining non-bonded interactions
NBONDS GROUP SWITCH CDIE VDW VSWI EXTEND GRAD QUAD -
CUTNB 22.0 CTOFNB 20.0 CTONNB 18.0 WMIN 1.5 EPS 1.0

!Calculating piston mass statistics
SCALar MASS STAT
CALCulate pmass = int (?STOT / 50.0)
CALCulate tmass = @pmass * 10

!Adaptive umbrella sampling for distance
rxncor: define set1 point sele atom PI 1 SOD end
rxncor: define set2 point sele atom NI 1 CLA end
rxncor: define rc distance set1 set2
rxncor: set rc
rxncor: trace rc unit 31
rxncor: UMBRella name rc KUMB @kumb DEL0 @del0 FORM 1
calc lowdel = @del0 - 1.5
calc hidel = @del0 + 1.5
rxncor: statistics name rc lowdelta @lowdel hidelta @hidel deldel 0.05 start 500

!Heating the system using thermostat
OPEN READ FORMatted UNIT 21 NAME @equ.res
OPEN WRITe FORMatted UNIT 20 NAME @umb.res
OPEN WRITe UNFOrmatted UNIT 22 NAME @umb.dcd

DYNA LEAP STARt -
NSTEP 150000 TIMESTep 0.002 NPRINT 1000 -
HOOVER REFT 300.0 TMASS @tmass
! PConst pmass 400.0 pgamma 20.0 tbath 300.0 PREFerence 1.0
IUNREAd -1 IUNWRIte 20 IUNCRD 22 NSAVC 50 -
IEQFRQ 0 IPRFRQ 1000 NTRFRQ 1000 INBFRQ -1 IMGFRQ -1 -
IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 1 ECHECK 2000

When I run using the above given input file, my output prints the following.

AVER DYN: Step Time TOTEner TOTKe ENERgy TEMPerature
AVER PROP: GRMS HFCTote HFCKe EHFCor VIRKe
AVER INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
AVER EXTERN: VDWaals ELEC HBONds ASP USER
AVER IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElec
AVER MMFP: GEO MDIP SSBP SHEL DROFfa
AVER PRESS: VIRE VIRI PRESSE PRESSI VOLUme
AVER UMBR: ADUMB RXNCor
---------- --------- --------- --------- --------- ---------
AVER> 1000 410.00000 -1728.20443 761.54885 -2489.75329 290.65515
AVER PROP> 19.18438 -1693.69510 865.06832 34.50934 870.21533
AVER INTERN> 322.74996 168.45926 0.00000 0.00000 0.00000
AVER EXTERN> 429.86435 -3599.40787 0.00000 0.00000 0.00000
AVER IMAGES> 0.00000 0.00000 0.00000 0.00000 -0.30749
AVER MMFP> 1.04264 0.00000 187.66405 0.00000 0.00000
AVER PRESS> -21.01478 -559.12877 0.00000 0.00000 0.00000
AVER UMBR> 0.00000 0.18182


How can the pressure be maintained at 0 atm?
There is something seriously wrong with my system and I need help to set up constant pressure and temperature for my system properly.



Edited by vidhya (09/07/18 06:24 AM)

Top
#37148 - 09/07/18 11:12 AM Re: Pressure unit in charmm [Re: Giselle]
rmv Online   content

Forum Member

Registered: 09/17/03
Posts: 8259
Loc: 39 03 48 N, 77 06 54 W
The HOOVER thermostat is part of the CPT code, and is probably not appropriate here.

I suspect the PRESS line in the output may not be valid; I believe it is zero because there is no defined unit cell volume.

Have you read the papers cited in the SSBP description of mmfp.doc?

I've never used SSBP myself, and that would be one of my next steps.
_________________________
Rick Venable
computational chemist


Top
#37155 - 09/18/18 02:42 AM Re: Pressure unit in charmm [Re: Giselle]
vidhya Offline
Forum Member

Registered: 08/09/18
Posts: 10
Loc: India
I read the paper cited in the mmfp.doc and they have used Langevin dynamics to control the temperature of the system.
I will try rerunning my simulations using LD and update.

Thanks
Vidhya

Top
Page 1 of 2 1 2 >

Moderator:  John Legato, lennart