You are not logged in. [Log In] CHARMM Development Project » Forums » User Discussion & Questions » Setup, I/O, and Basic questions » Pressure unit in charmm Register User    Forum List        Calendar     Active Topics    Search     FAQ
 Page 1 of 2 1 2 >
 Topic Options
 #25199 - 08/27/10 03:23 AM Pressure unit in charmm Giselle 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 360I 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 VOLUmeLAVE PRESS> 478.01444 -1155.66231 -1112.90548 23.11227 29451.39517The 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^21 N = 1.4388 kcal/mol-A1 m^2 = 1E20 A^21 atm = 101325*1.4388*1E-20 kcal/mol-A^3 =1.4579E-15 kcal/mol-A^3Are 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 lennart Forum Member Registered: 09/25/03 Posts: 4761 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 NilssonKarolinska InstitutetStockholm, Sweden Top
 #25201 - 08/27/10 05:35 AM Re: Pressure unit in charmm Giselle 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 lennart Forum Member Registered: 09/25/03 Posts: 4761 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 NilssonKarolinska InstitutetStockholm, Sweden Top
 #25204 - 08/27/10 06:21 AM Re: Pressure unit in charmm Giselle 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 rmv Forum Member Registered: 09/17/03 Posts: 8419 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 Venablecomputational chemist Top
 #25217 - 08/28/10 01:12 AM Re: Pressure unit in charmm [Re: rmv] Kenno 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 vidhya Forum Member Registered: 08/09/18 Posts: 13 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 potentialMMFPSSBP SELE TYPE OH2 ENDGEO sphere RCM - xref 0.0 yref 0.0 zref 0.0 - force 10.0 droff 0.0 - select segid PI .or. segid NI endEND!Defining non-bonded interactionsNBONDS 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 statisticsSCALar MASS STATCALCulate pmass = int (?STOT / 50.0)CALCulate tmass = @pmass * 10!Adaptive umbrella sampling for distancerxncor: define set1 point sele atom PI 1 SOD endrxncor: define set2 point sele atom NI 1 CLA endrxncor: define rc distance set1 set2rxncor: set rcrxncor: trace rc unit 31rxncor: UMBRella name rc KUMB @kumb DEL0 @del0 FORM 1calc lowdel = @del0 - 1.5calc hidel = @del0 + 1.5rxncor: statistics name rc lowdelta @lowdel hidelta @hidel deldel 0.05 start 500!Heating the system using thermostatOPEN READ FORMatted UNIT 21 NAME @equ.resOPEN WRITe FORMatted UNIT 20 NAME @umb.resOPEN WRITe UNFOrmatted UNIT 22 NAME @umb.dcdDYNA 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 2000When I run using the above given input file, my output prints the following.AVER DYN: Step Time TOTEner TOTKe ENERgy TEMPeratureAVER PROP: GRMS HFCTote HFCKe EHFCor VIRKeAVER INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopersAVER EXTERN: VDWaals ELEC HBONds ASP USERAVER IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElecAVER MMFP: GEO MDIP SSBP SHEL DROFfaAVER PRESS: VIRE VIRI PRESSE PRESSI VOLUmeAVER UMBR: ADUMB RXNCor ---------- --------- --------- --------- --------- ---------AVER> 1000 410.00000 -1728.20443 761.54885 -2489.75329 290.65515AVER PROP> 19.18438 -1693.69510 865.06832 34.50934 870.21533AVER INTERN> 322.74996 168.45926 0.00000 0.00000 0.00000AVER EXTERN> 429.86435 -3599.40787 0.00000 0.00000 0.00000AVER IMAGES> 0.00000 0.00000 0.00000 0.00000 -0.30749AVER MMFP> 1.04264 0.00000 187.66405 0.00000 0.00000AVER PRESS> -21.01478 -559.12877 0.00000 0.00000 0.00000AVER UMBR> 0.00000 0.18182How 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 rmv Forum Member Registered: 09/17/03 Posts: 8419 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 Venablecomputational chemist Top
 #37155 - 09/18/18 02:42 AM Re: Pressure unit in charmm vidhya Forum Member Registered: 08/09/18 Posts: 13 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.ThanksVidhya Top
 #37163 - 09/26/18 01:52 AM Re: Pressure unit in charmm vidhya Forum Member Registered: 08/09/18 Posts: 13 Loc: India HiI performed the MD run using Langevin dynamics as given in the CHARMM tutorial.I am able to run the simulation without any issues in temperature control. My only concern is the pressure.The log file shows that the pressure is O atm.This means the simulation is running as such it is in vacuum.I have given the scalar fbeta to be 5 ps-1 as per the reference paper given in mmfp.doc.From what I read, LD is mostly used for implicit solvent simulations. And we do not define the pressure but still LD says it runs in NPT ensemble.Here is a sample of my output:DYNA DYN: Step Time TOTEner TOTKe ENERgy TEMPeratureDYNA PROP: GRMS HFCTote HFCKe EHFCor VIRKeDYNA INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopersDYNA EXTERN: VDWaals ELEC HBONds ASP USERDYNA IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElecDYNA MMFP: GEO MDIP SSBP SHEL DROFfaDYNA PRESS: VIRE VIRI PRESSE PRESSI VOLUmeDYNA UMBR: ADUMB RXNCor GAMUS ---------- --------- --------- --------- --------- ---------DYNA> 2000 4.00000 -1630.14038 808.23955 -2438.37994 307.77499DYNA PROP> 21.74849 -1585.83007 937.07987 44.31031 1170.23642DYNA INTERN> 394.96040 168.31806 0.00000 0.00000 0.00000DYNA EXTERN> 425.32588 -3617.00872 0.00000 0.00000 0.00000DYNA IMAGES> 0.00000 0.00000 0.00000 0.00000 0.62431DYNA MMFP> 0.20326 0.00000 189.03856 0.00000 0.00000DYNA PRESS> -18.12466 -762.03295 0.00000 0.00000 0.00000DYNA UMBR> 0.00000 0.15832 0.00000 ---------- --------- --------- --------- --------- ---------My pressure is 0 throughout the simulation.How can I set the pressure constant to be around 1 atm?Input file:!Defining spherical solvent boundary potentialMMFPSSBP SELE TYPE OH2 ENDGEO sphere RCM - xref 0.0 yref 0.0 zref 0.0 - force 10.0 droff 0.0 - select segid PI .or. segid NI endENDCONS HARM BESTFIT MASS FORCE 10.0 SELE SEGId PI END NOROTCONS HARM BESTFIT MASS FORCE 10.0 SELE SEGId NI END NOROT!Initialized Langevin dynamicsSCALAR FBETA SET 5.0 SELEct TYPE OH2 .AND. SEGId W END!Defining non-bonded interactionsNBONDS GROUP SWITCH CDIE VDW VSWI EXTEND GRAD QUAD -CUTNB 22.0 CTOFNB 20.0 CTONNB 18.0 WMIN 1.5 EPS 1.0!Adaptive umbrella sampling for distancerxncor: define set1 point sele atom PI 1 SOD endrxncor: define set2 point sele atom NI 1 CLA endrxncor: define rc distance set1 set2rxncor: set rcrxncor: trace rc unit 31rxncor: UMBRella name rc KUMB @kumb DEL0 @del0 FORM 1calc lowdel = @del0 - 1.5calc hidel = @del0 + 1.5rxncor: statistics name rc lowdelta @lowdel hidelta @hidel deldel 0.05 start 500!Heating the system using thermostatOPEN READ FORMatted UNIT 21 NAME @equ.resOPEN WRITe FORMatted UNIT 20 NAME @umb.resOPEN WRITe UNFOrmatted UNIT 22 NAME @umb.dcdDYNA LEAP LANGEVIN STARt - NSTEP 150000 TIMESTep 0.002 NPRINT 1000 - FIRSTT 300.0 FINALT 300.0 TBATH 300.0 IHTFRQ 0 TWINDH 5.0 TWINDL -5.0 - IUNREAd -1 IUNWRIte 20 IUNCRD 22 NSAVC 50 - IEQFRQ 0 IPRFRQ 1000 NTRFRQ 1000 INBFRQ -1 IMGFRQ -1 - IASORS 1 IASVEL 1 ISCVEL 0 ICHECW 1 ECHECK 2000 Top
 #37164 - 09/26/18 10:29 AM Re: Pressure unit in charmm rmv Forum Member Registered: 09/17/03 Posts: 8419 Loc: 39 03 48 N, 77 06 54 W The output also shows a volume of zero; the pressure cannot be computed w/o a defined volume. The simulation pressure codes in CHARMM require the use of periodic boundary conditions.You may be able to get a crude estimate of the pressure in an ad hoc way, before or after a simulation, using the PRESS INST command and supplying a volume and temperature. Estimating the volume of your droplet may be challenging. You may wish to research the topic of computing pressure in a droplet.You could also try contacting the SSBP code authors. _________________________ Rick Venablecomputational chemist Top
 Page 1 of 2 1 2 >

 Hop to: User Discussion & Questions ------   Setup, I/O, and Basic questions   Energy terms, Constraints, Restraints, and Solvation   Parameter Set Discussion   Molecular Dynamics   Minimization, Normal modes, Monte Carlo,...   Installation and Testing   Performance and Benchmarking   Other User Discussion and Questions   General Chemistry DiscussionsCHARMM Interfaces ------   QM/MM Discussion and Questions   CHARMMing   CHARMM-GUI   MMTSB   AccelrysCHARMM Community ------   Script Archive   Bug Reports & Fixes   Meetings & Events   Position Openings   CHARMM Course   Suggestions and Requests   Forum News and Announcements
Moderator:  John Legato, lennart