CHARMM Development Project
Posted By: Allen_123 Post Analysis of openMM trajctory - 01/14/21 05:01 AM
Hello All,

Is there any way Charmm can read pdb and dcd produced by openMM? The code that Charmm reads pdb and dcd cannot open files of OpenMM. I know there are software like mdtraj for post analysis of trajectory, but only Charmm has an option to measure nonbond interactions between different parts of molecule.

OPEN UNIT 1 FORM READ NAME step3_pbcsetup.psf
READ PSF CARD UNIT 1
CLOSe UNIT 1
OPEN UNIT 1 READ CARD NAME PFF_8cell_last.pdb
READ COOR pdb UNIT 1
CLOSE UNIT 1

open unit 21 read unform name PFF_8cell_10ns.dcd

OPNLGU> Unit 1 cannot be opened as PFF_8CELL_LAST.PDB

***** LEVEL 0 WARNING FROM *****
***** "OPEN" not possible.
Posted By: Allen_123 Re: Post Analysis of openMM trajctory - 01/14/21 07:37 AM
I figured out this error by changing the file name to pff_8cell_10ns.dcd
Posted By: lennart Re: Post Analysis of openMM trajctory - 01/14/21 08:39 AM
The following should also work (note that there is no need to have OPEN and CLOSE commands with READ and WRITE commands):
read psf card name step3_pbcsetup.psf
read coor pdb name "PFF_8cell_last.pdb"
open unit 21 read unform name "PFF_8cell_10ns.dcd"

Using double quotes preserves case in filenames when you read files.
Posted By: Allen_123 Re: Post Analysis of openMM trajctory - 01/14/21 11:20 AM
Thank you for your fast respond.

Another question is how to extract energy from interaction command? For instance, in the following commands, I want to write vdw and elec energy to a file.

open read unit 21 file name pff_8cell_10ns.dcd
traj query unit 21
traj iread 21 nunit 1
set x 1
label loop
traj read
energy
update
SKIP ALL EXCL vdw Elec
interaction select resname phe .and. type C* end sele resname tip3 .and. type OH2 end
Posted By: rmv Re: Post Analysis of openMM trajctory - 01/14/21 04:27 PM
The energy doc file has a detailed list of energy terms available via ? substitutions, including ?ELEC and ?VDW, which can be used to write those values to a file:

open read unit 21 file name pff_8cell_10ns.dcd
OPEN UNIT 2 WRITE CARD NAME INTE.TXT
ECHU 2 ! see miscom doc file for echo usage
traj query unit 21
SET COUNT = ?NFILE
traj iread 21 nunit 1
set x 1
label loop
traj read
energy
update
SKIP ALL EXCL vdw Elec
interaction select resname phe .and. type C* end sele resname tip3 .and. type OH2 end
ECHO @X ?VDW ?ELEC
INCR X BY 1
IF X LE @COUNT GOTO LOOP


A more complete list of ? substitutions can be found in the subst doc file. Also, there is no need to read the PDB file for this, the PSF is sufficient for reading DCD files.
© CHARMM forums