Previous Thread
Next Thread
Print Thread
Joined: Mar 2005
Posts: 88
R
rossi Offline OP
Forum Member
OP Offline
Forum Member
R
Joined: Mar 2005
Posts: 88
Good Morning,

I am simulating a system (description below) using information derived from J.Chem. Phys. 103, 10252 (1995) and J.Chem. Phys. 103, 10267 (1995).

The system is a water/heptane system, constructed with the heptane placed in the center of the box and the water above and below. Each cube is 54 A on a side.

NVT Simulations
=============

! set up nonbond parameters -- same as for heating
nbond inbfrq -1 imgfrq -1 atom vatom cdie eps 1.0 -
elec ewald pmew fftx 64 ffty 64 fftz 64 kappa .34 spline order 6 -
vdw vswitch cutnb 16.0 cutim 16.0 ctofnb 12. ctonnb 10. wmin 1.0

and

! NVT Dynamics Command
dyna leap cpt start - ! restart the run from our heating restart, use CPT dynamics
timestep 0.001 nstep 100000 nprint 1000 - ! do 100.0 ps in 1 fs time step
iunrea -1 iunwri 41 iuncrd 51 isvfrq 1000 nsavc 1000 - ! units for reading/writing restarts and coordinate trajectory
hoover reft 300.0 tmass @tmass - ! Hoover thermostat, constant temp - 300 K, tmass as above
pcons pref 1.0 pmass 0.0 pgamma 20 - ! NVT simulation, pmass 0.5 (need pmass~0), collision frequency 0 ?????
ihtfrq 0 ieqfrq 0 ntrfrq 1000 - ! temperature under Hoover's control, cancel rotation/translation every 1000 steps
iasors 1 iasvel 1 iscvel 0 ichecw 0 twindh 5.0 twindl -5.0 - ! since we're restarting and ihtfrq = ieqfrq = 0, these don't really matter
echeck 100. ! check energy diff as before

A surface tension value is printed out every 1000 ts. Although I haven't tried it, I would expect an NVT simulation on either the separate water or heptane boxes would yield a surface tension value ~ 0 because these systems are isotropic and have images.

My question is related to the system with two water boxes and the heptane box. Since there are two interfacial regions, should I divide the value for surface tension provided by the NVT simulation by 2? My reason for asking is that for the NPAT simulation for the octane/water system described in J.Chem. Phys. 103, 10252 (1995) where the surface tension, gamma = 1/2 h_z (p_zz - (p_xx + p_yy)/2) . Here the 1/2 is required
because of the two interfacial regions.

Could someone please explain to me the algorithm used in computing the surface tension for an NVT simulation. I would really like that.

NPAT Simulations
==============

I just started those simulations and a similar set of commands are used, but pmass 0.0 is changed tp pxx = pyy = 0 and pzz = zzzz

Any comments about these simulations would be very helpful. Thanks so much for your help with this.

Regards,

Angelo


P.S. How does one carry out NPxH type simulations discussed in J.Chem. Phys. 103, 10252 (1995) using charmm? I am sure it is straight forward, but I am clueless about it. This is only for my own sense of understanding the papers.

Joined: Sep 2003
Posts: 8,605
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
The implementation assumes two interfaces, and reports the tension per interface. (It is more correctly called the interfacial tension, where surface is reserved for liquid:air or liquid:solid systems.)

NPH indicates a barostat, but no thermostat; it is rarely used.

Due to the periodic boundary conditions, I believe your system has one heptane phase and one water phase; the description is a bit vague.


Rick Venable
computational chemist

Joined: Mar 2005
Posts: 88
R
rossi Offline OP
Forum Member
OP Offline
Forum Member
R
Joined: Mar 2005
Posts: 88
Rick,

There is a box of heptane molecules (50Ax50Ax50A), equilibrated for about 5 ns, in the center of the unit cell

On the TOP of that box of heptane molecules is a box of water molecules (50Ax50Ax50A), equilibrated for about 5 ns ---> this constitutes ONE heptane/water surface for interfacial surface tension.

On theBOTTOM of that box of heptane molecules is a box of water molecules (50Ax50Ax50A), equilibrated for about 5 ns ---> this constitutes A SECOND heptane/water surface for interfacial surface tension.

It is the same construct used for the water/octanol system used in the papers mentioned:

----------------
| H2O |
-------------- ----------> interfacial interaction and surface tension
|HEPTANE|
-------------- ----------> interfacial interaction and surface tension
| H2O |
----------------

Joined: Sep 2003
Posts: 8,605
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
The original water:octane simulations in the paper had a simulation box that was was twice as long in z as it was in x and y, i.e. the initial volumes were roughly the same for both fluids. Your system is 3 times longer in z, and has twice the volume of water as heptane, so it's not quite the same, but the extra water should not have much impact on the results.


Rick Venable
computational chemist

Joined: Mar 2005
Posts: 88
R
rossi Offline OP
Forum Member
OP Offline
Forum Member
R
Joined: Mar 2005
Posts: 88
OK. Thank you, Rick.

Joined: Mar 2005
Posts: 88
R
rossi Offline OP
Forum Member
OP Offline
Forum Member
R
Joined: Mar 2005
Posts: 88
Rick,

I would like to follow up on determining the surface tension using the NVT simulation input I descrbed previously in this post. Here is the result for one of the iterations:

DYNA DYN: Step Time TOTEner TOTKe ENERgy TEMPerature
DYNA PROP: GRMS HFCTote HFCKe EHFCor VIRKe
DYNA INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
DYNA EXTERN: VDWaals ELEC HBONds ASP USER
DYNA IMAGES: IMNBvdw IMELec IMHBnd RXNField EXTElec
DYNA EWALD: EWKSum EWSElf EWEXcl EWQCor EWUTil
DYNA PRESS: VIRE VIRI PRESSE PRESSI VOLUme
DYNA XTLE: XTLTe SURFtension XTLPe XTLtemp
---------- --------- --------- --------- --------- ---------
DYNA> 7000 73.00000-305799.21877 32723.67828-338522.89706 297.89708
DYNA PROP> 14.46592-305764.41471 32828.10517 34.80406 37909.45962
DYNA INTERN> 755.49207 4399.76220 1327.21100 2968.63336 0.00000
DYNA EXTERN> 10147.88739-172730.47223 0.00000 0.00000 0.00000
DYNA IMAGES> -456.08333 -11336.55832 0.00000 0.00000 0.00000
DYNA EWALD> 4146.03374-923597.89259 745853.08965 0.00000 0.00000
DYNA PRESS> -4265.23656 -21007.73652 520.69670 107.14488 561670.27012
DYNA XTLE> -361809.62597 282.82939 -56010.69932 0.00000

The surface tension output is 282.82939. Values for water/alkane interfacial tension for heptane/water should be ~46 dyne/cm

What are the units? Are they dyne/cm?

Well, I assumed the output values are not converged.

But, a continuing simulation seems to converge to around a value of 220.

What are the units, please? How is surface tension computed in NVT charmm calculations?

One final disparate comment, and I will open a separate post if you think I need to.
Parallel performance of the free version of charmm is terrible! I understand there is a more efficient version out there. But, assuming all my years in the Computer Science Department working on HPC applications is worth something, I believe that the vertical/lateral parallel scaling of the free version is really poor. The best performance seems to be with a single node multicore mode of 24 cores. Even moving up to a single node with 36 cores really doesn't help at all. And trying more than one node is impossible, i.e. really doesn't work. Something has happened, and I can't understand it. In the past, all version of charmm scaled firly well. Not any more. It really hampers my own research output. There is a definition of a supercomputer we used to use at the IBM Research Division at Watson: a supercomputer is one that increases the publication of quality research articles. charmm free version doesn't meet that definition. Gosh, I rreally love charmm. It's my goto program. What's happened to it?

Sorry for the rant!

Regards,

Angelo

Joined: Sep 2003
Posts: 8,605
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
Dyne/cm

Simulations of many ns are needed to get converged average values for surface tension, as the pressure fluctuations in MD systems are substantial. Snapshots values are not expected to be very close to the target. Run for 100 ns, and look at the time series. Maybe read some of our papers that report interfacial tension for lipid systems.

The non-free (still only $600) version includes the fast MD engine (domdec), which scales comparable to NAMD when identical conditions are used, esp. timestep and PME grid update frequency. NAMD is faster with thousands of cores, but the thermostats and barostats are not as rigorous as those in CHARMM; they made that sacrifice for speed.

I use OpenMM (not in CHARMM) for many of my simulations now; it's fast GPU code that gives ensembles that match those from CHARMM, esp. with the new Nose-Hoover thermostat. The exception is simulations that need more rigorous pressure calculations, such as interfacial tension, viscosity, etc. The virial-based pressure from every integration step used in CHARMM is the only suitable approach; most fast MD programs use a Monte Carlo barostat, which cannot be used for those calculations.

Questions about features in free and non-free should be directed to Prof. Martin Karplus.


Rick Venable
computational chemist

Joined: Mar 2005
Posts: 88
R
rossi Offline OP
Forum Member
OP Offline
Forum Member
R
Joined: Mar 2005
Posts: 88
Rick, thanks so much for your patience. In addition, part of my frustration is that my most of my jobs are queued on the UConn cluster. In any case I will look into OpenMM.
I should mention, that charmm is really great for setting up complicated input files for projects. And I really like charmm for analysis. It simply cannot be beat for that. The requirements are usually a psf file, a crd or pdb file, and dcd files, and you can then do almost any analysis possible. The script library is excellent! Thanks so much again, Angelo

Joined: Sep 2003
Posts: 4,850
Likes: 10
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 4,850
Likes: 10
Note that you can use OpenMM from CHARMM (even with the free version). Check doc/openmm.info and the testcases test/c*test/omm*


Lennart Nilsson
Karolinska Institutet
Stockholm, Sweden
Joined: Sep 2003
Posts: 8,605
Likes: 24
rmv Online Content
Forum Member
Online Content
Forum Member
Joined: Sep 2003
Posts: 8,605
Likes: 24
For my work, I've needed some of the newer OpenMM features not yet available in the CHARMM interface, esp. the new, more efficient Nose-Hoover thermostat with multiple heat baths (for Drude force fields), as was as the LJ-PME method for long range VDW. Fortunately, one of our former post-docs put together a python package (install required) that provides a single, simple python script that has all of the same information you would put in a CHARMM input script. I have had to improve my knowledge of python, which has been worth the effort.

I still use CHARMM for most analyses, using the DCD format files from OpenMM, and the same PSF I used to build the initial model with CHARMM.


Rick Venable
computational chemist


Moderated by  BRBrooks, lennart, rmv 

Link Copied to Clipboard
Powered by UBB.threads™ PHP Forum Software 7.7.5
(Release build 20201027)
Responsive Width:

PHP: 7.3.31-1~deb10u1 Page Time: 0.015s Queries: 34 (0.010s) Memory: 0.7836 MB (Peak: 0.8686 MB) Data Comp: Off Server Time: 2022-09-28 01:12:23 UTC
Valid HTML 5 and Valid CSS