From: Prathit Chatterjee (pc20apr_at_REMOVE_yahoo.co.in)
Date: Fri Nov 20 2020 - 04:49:42 CST

Dear All,

Here is again a related problem to the previous chain of emails.

While I can load a netcdf trajectory generated with new AMBER versions, I cannot write a modified netcdf AMBER coordinate using VMD, instead I can only write as a crd/crdbox format, please see below kindly:

vmd > package require pbctools2.8vmd > mol new /homes/epsilon/users/diem/project_tau43_Ab42/tau43_monomer/traj_10/production/PHF43.top type parm7Info) Using plugin parm7 for structure file /homes/epsilon/users/diem/project_tau43_Ab42/tau43_monomer/traj_10/production/PHF43.topInfo) Analyzing structure ...Info)    Atoms: 144684Info)    Bonds: 144686Info)    Angles: 0  Dihedrals: 0  Impropers: 0  Cross-terms: 0Info)    Bondtypes: 0  Angletypes: 0  Dihedraltypes: 0  Impropertypes: 0Info)    Residues: 36056Info)    Waters: 36009Info)    Segments: 1Info)    Fragments: 36014   Protein: 1   Nucleic: 00vmd > mol addfile /homes/epsilon/users/diem/project_tau43_Ab42/tau43_monomer/traj_10/production/crd/PHF43_310K_10_md020.crd waitfor all type netcdfnetcdfplugin) conventions: 'AMBER'netcdfplugin) trajectory follows AMBER conventions version '1.0'netcdfplugin) AMBER: program 'pmemd'netcdfplugin) AMBER: program version '16.0'netcdfplugin) AMBER: title 'default_name'netcdfplugin) AMBER: application 'AMBER'netcdfplugin) AMBER: spatial dimension: 3netcdfplugin) AMBER: atom dimension: 144684netcdfplugin) AMBER: frame dimension: 200netcdfplugin) AMBER: coordinates units: 'angstrom'netcdfplugin) AMBER: no coordinates scalefactor attribute, 1.0 assumednetcdfplugin) AMBER: coordinates scalefactor: 1.000000netcdfplugin) AMBER trajectory contains periodic cell informationnetcdfplugin) AMBER: cell lengths units: 'angstrom'netcdfplugin) AMBER: no cell lengths scalefactor attribute, 1.0 assumednetcdfplugin) AMBER: cell lengths scalefactor: 1.000000netcdfplugin) AMBER: cell angles units: 'degree'netcdfplugin) AMBER: no cell angles scalefactor attribute, 1.0 assumednetcdfplugin) AMBER: cell angles scalefactor: 1.000000Info) Using plugin netcdf for coordinates from file /homes/epsilon/users/diem/project_tau43_Ab42/tau43_monomer/traj_10/production/crd/PHF43_310K_10_md020.crdInfo) Finished with coordinate file /homes/epsilon/users/diem/project_tau43_Ab42/tau43_monomer/traj_10/production/crd/PHF43_310K_10_md020.crd.0vmd > set all [atomselect top all]atomselect0vmd > pbc wrap -centersel "protein" -center com -compound residue allerror: pbcwrap: unknown option: allvmd > pbc wrap -centersel "protein" -center com -compound residue -allInfo) 0.5% complete (frame 0)Info) 18.5% complete (frame 36)Info) 38.0% complete (frame 75)Info) 57.5% complete (frame 114)Info) 77.0% complete (frame 153)Info) 96.0% complete (frame 191)Info) 100.0% complete (frame 199)vmd > animate write netcdf PHF43_310K_10_md020.crd beg 0 end 199 waitfor all sel $allERROR) Unable to open file PHF43_310K_10_md020.crd of type netcdf for writing frames.vmd > animate write crdbox wrap_310K_10_md020.crd beg 0 end 199 waitfor all sel $allInfo) Opened coordinate file wrap_310K_10_md020.crd for writing.Info) Coordinate I/O rate 5.99 frames/sec, 9 MB/sec, 33.4 secInfo) Finished with coordinate file wrap_310K_10_md020.crd.200
In order to avoid losing resolution by not writing coordinates in ASCII format, any suggestions to write coordinate files from the original one in the same NETCDF format, will be deeply appreciated.
Thanking you in advance,
Prathit Chatterjee

    On Monday, 9 November, 2020, 06:22:17 pm GMT+9, Peter Freddolino <petefred_at_umich.edu> wrote:
 
 I think one useful step would be to actually calculate the magnitudes of the differences in coordinates between the original trajectory and the one read/written by VMD. The energy differences that you've noted here are tiny, on the order of less than one part per thousandth, and the energies are going to be very sensitive to small changes in coordinates.
I also suspect that if you actually inspect the files before/after running them through your procedure in vmd, you will find that they are in fact not the same format. Since you're using netcdf files in amber, you have a binary container format when you read them into VMD, but when you write from VMD in crdbox format, you're converting to an ascii format with different precision. Just because the file extension is .crd doesn't actually mean that the format is the same; you can see this easily enough by looking at them in any text viewer or editor. This seems likely to me to be the issue. If you want to avoid any changes you'd be best off comparing apples to apples, and having either netcdf-format or crdbox-format files both before and after vmd sees the coordinates.Best,Peter

On Sun, Nov 8, 2020 at 10:56 PM Prathit Chatterjee <pc20apr_at_yahoo.co.in> wrote:

Dear Prof. Freddolino,

Sorry for the delayed response.

The precision of selected coordinate numbers to calculate the energy is the same in both the original AMBER coordinates and protein-only generated trajectory (in VMD), and the format is .crd, calculated with cpptraj, please see below:

This is for the original amber coordinates:
parm prot.top
trajin /homes/epsilon/users/diem/project_t-ab/traj_1/production/crd/T_310K_1_md001.crd 1 10

#centering
strip :WAT
strip :Na+
strip :Cl-

center :1-43 mass origin
image origin center familiar

pairwise IE :1-43 out pairwise_IE.out cuteelec 0.0001 cutevdw 0.0001 vmapout IE_vmap.out emapout IE_emap.out avgout IE_avg.out
--------------

This is for the protein only generated coordinates:parm prot.top
trajin /homes/epsilon/users/pchatterjee/vmd_anal_prep/diem_1/T_310K_prot_diem_1_100ns.crd 1 10

pairwise IE :1-43 out pairwise_IE.out cuteelec 0.0001 cutevdw 0.0001 vmapout IE_vmap.out emapout IE_emap.out avgout IE_avg.out
-------------

Also, I generated the protein-only amber coordinates as follows in VMD (a prototypical example provided below):
        mol new /homes/epsilon/users/diem/project_t-ab/traj_10/production/whole.top type parm7
        mol addfile /homes/epsilon/users/diem/project_t-ab/traj_10/production/crd/T_310K_10_md$n.crd type netcdf  waitfor all

        set frames [molinfo top get numframes]
        set nf [expr {[molinfo top get numframes]-1}]

        set a [atomselect top protein]
        animate write crdbox diem_10/PHF43_310K_prot_$f.crd beg 0 end $nf waitfor all sel $a

I tried to write the protein only coordinates in several ways until reaching the above mentioned procedure, loading the original coordinates as netcdf, and writing the final protein only coordinates as crdbox.
Likewise, the coordinates of the original and protein-only coordinates could be visualized properly without noticing any apparent unnatural bonds.
Best Regards,
Prathit On Sunday, 8 November, 2020, 01:35:42 pm GMT+9, Peter Freddolino <petefred_at_umich.edu> wrote:
 
 Dear Prathit,We really need more information on how you did these calculations. How were the energies calculated? And how exactly did you generate the trajectory in vmd? With such small differences it could be something as simple as a difference in the numerical precision at which coordinates were stored, depending on exactly what formats you used.Best,
Peter

On Sat, Nov 7, 2020 at 10:02 AM Prathit Chatterjee <pc20apr_at_remove_yahoo.co.in> wrote:

 Dear VMD experts,
I had compiled a protein only amber trajectory of an explicitly solvated simulation with VMD, which I was able to visualize properly and run certain analyses with it in VMD as well.
But, for cross-checking, I tried to calculate the pairwise interaction of the original AMBER coordinates, and  the protein only generated trajectory (with VMD).
The results, although very similar, are showing certain discrepancy as follows:
[pchatterjee_at_node51 pairwise_IE]$ tail -f test*/pairwise_IE.out
==> test1/pairwise_IE.out <== [THIS IS THE ORIGINAL AMBER COORDINATE]
       1    -186.3642   -2973.0974
       2    -189.9812   -2966.1259
       3    -199.5370   -2957.7041
       4    -196.3846   -2983.7931
       5    -197.8264   -2981.1530
       6    -184.0711   -2971.6430
       7    -194.3659   -2970.1477
       8    -189.0441   -2956.9308
       9    -197.0864   -2948.2561
      10    -196.4939   -3010.1367

==> test2/pairwise_IE.out <==[THIS IS THE PROTEIN ONLY GENERATED AMBER COORIDNATE IN VMD]
       1    -186.3943   -2972.9815
       2    -189.9798   -2966.1573
       3    -199.5455   -2957.6502
       4    -196.3961   -2983.7881
       5    -197.8009   -2981.1609
       6    -184.0410   -2971.7176
       7    -194.3681   -2970.0830
       8    -189.0185   -2956.9181
       9    -197.1040   -2948.2302
      10    -196.4685   -3010.0947

Although the global picture will be the same with such minor discrepancy, I just want to know for the sake of knowing what the underlying reason could be.
Any comments/suggestions will be extremely helpful.
Thank you in advance,Best Regards,Prathit