Klauda Lab Simulation Wiki

First of all, please check the dcdfreq and xstfreq in your dyn-2.inp file, if the dcdfreq is equal to xstfreq (which is the case almost all the time), then do A and skip B; if the dcdfreq is not equal to xstfreq, then skip A and do B only. The A and B session both generate the *-cut.dat. (Since the coord of each lipid is based on each frame, so we want to have the number of x-dimension data point agree with number of frame. You may know this better after you finish this rdf analysis.)

A. Update the trim.scr

# ./trim.csh (first dcd) (last dcd) (path of xst file)
 
# UPDATE the three variables mentioned above
set f = 10
set l = 11
set p = /lustre/jbklauda/lipid/analysis-test/gl-chl1

Note: The p is the charmm-gui path that has the dyn*.xst.gz files and f and l are the first and last dyn*.dcd file.

Run trim.scr by command:

./trim.scr

This will get the x-dimension data for each frame. (If you do A, you don’t need the extract folder, please ignore it.) The dyn*-cut.dat will be generated in rdf folder after run.

B. Make sure you have copied the extract folder and then go to extract folder.

Update the extract.scr

# ./extract.csh (first dcd) (last dcd) (path of psf file)
 
# UPDATE the three variables mentioned above
set f = 1
set l = 75
set p = /lustre/aou1/soybean/hypocotyl-b1

Note: The p is the charmm-gui path that has the step5_assembly.psf file and f and l are the first and last dyn*.dcd file.

Update the path of topology files in gen_crd.inp and extract.inp if needed.

Run extract.scr by command:

extract.scr

This does two steps, 1) extract the crd file (The extract crd is used as reference crd file instead of step5_assembly.crd, just to make the image of lipid agree in all dcd.) 2) x-dimension data for each frame is extracted. The dyn*-cut.dat will be generated in rdf folder after run.

1. Copy the toppar.str and toppar subdirectory from the CHARMM-GUI folder (CHARMM input files) to the 2d-rdf directory. Also, update the psf and crd file in psfcrd.str. This must be an absolute path so do not use relative paths for the psf/crd files.

2. In 2d-rdf folder, make a directory called “inp” if you don’t have it (to store the input file that will be read into the *.exe program when calculating the RDF between lipids).

3. Update the size and data point information in gen-inp.csh if necessary. Obtain the x dimension value of the last value from dyn.xsc (in your dyn folder), first number is about 2/3 of x-dim in simulation box e.g. if the value is 72.1711, then 72*2/3=48≈50, then 50 is used as the first number.

# first number is ~2/3 of x-dim in simulation box 
# second number is number of points in g(r)
# ~0.1 incremets across (r) so funct looks smooth
  echo "50 500" >> inp/inp$k

In the example, assure two third of your system size is about 50 Å (the first value above), then 50/0.1=500, so you use 500 points (the second value above), if it is 60 Å, then use 600 points, etc.

4. Update the dyn file range in gen-inp.scr and then run command ./gen-inp.scr to generate input file for RDF calculation.

You will have many pairs of heagroup that you want to calculate based on your system. Please work on the folders with same lipid headgroup first, after they are all finished, then work on the different (mixed) lipids because we need to find the coordinate of each lipid first, and then the mixed lipids pair just using the coordinated of lipid head group already calculated. In the example, after finish the calculation in same lipid pair, e.g. pc,pe,pi,sterols, then we start to work on pc-pe, pc-pi,pc-sterol, etc.

The scripts for same lipid pair for glycerol (please use eg_same-glycerol as script template) and sterol lipids (eg_same-sterol), and also the mixed lipids pairs are slightly different (eg_diff). So for each lipid headgroup pair, please copy the correct folder/ scripts as template and rename the folder based on the headgroup names in your system (e.g. use following command to copy the glycerol script folder for your pc lipid group: scp -r eg_same-glycerol pc ), and then do the following steps:

5. In the each of same (single) lipid folder (eg. pc,pe,pi,sterol etc.), Update the residue names and number of types of each headgroup in coor-*.inp

set rn1 = PLPE                                           ! Update the resname
set rn2 = LLPE                                           ! Update the resname
set rn3 = L2PE                                           ! Update the resname
 
set Ntyp = 3   ! number of types of lipids, which should be consistent as resnames types above

5.1 Update the job name in coord.csh and rdf.csh (This just to distinguish to make debugging easier if error occurs.) Also, copy the toppar.str and toppar subdirectory into the respective eg* directories.

5.2 Updating coord.scr by update the range and path of DCD files to be read. If needed, update the charmm version used in coor.csh.

5.3 Run command ./coord.scr

After run coor.scr, the x,y, and z coordinate of specified atom (P or C2 for glycerol lipids (selecting baesd on need), O3 for sterols, etc.) of each lipid in each frame are in xyz-*.dat.

5.4 Update the dcd range in rdf.scr as you will compute the RDF only for the equilibrated portion of the trajectory rdf.csh (first dcd) (last dcd) (first equilibrated dcd). The (first equilibrated dcd) can be the same as the (first dcd), but if you wan the rdfs for all dcds but only do averages for a equilibrated section then (first dcd) would be 1 and (irst equilibrated dcd) would be the first dcd in the equilbrated section. Once the values are entered run ./rdf.scr.

6. After the same headgroup calculation is finished, then we can start to work for the different (mixed) lipid folder(eg. pc-pe,pc-pi,pc-sterol etc.),

6.1 Update the lipid head group (folder) name in gen-in.scr, edit the range of r as step 3 of Part II, and then run gen-in.scr
6.2 Update the dcd range in rdf.scr as you will compute the RDF only for the equilibrated portion of the trajectory rdf.csh (first dcd) (last dcd) (first equilibrated dcd). And then run rdf.scr.

Note: You don’t have coor.csh nor coor.scr since the all the coordinate file that we need are already calculated in single lipid folder. All we need to do is to generate the new input files that select the specified coordinate files (in order to calculate 2d rdf in the next step). The gen-inp.csh for same lipid (one directory out) which select only one xyz-p*.dat file, which is different with gen-inp.csh in mixed lipid folder, which select two xyz-p*.dat files of different lipids.

7. If you are interested in the block 2d rdf instead of each dcd, update the dcd range in each block, and run block.scr, the block values will be shown in avg-rdf-p.dat. The number of block can be selected based on your interest. (e.g. 5 to 10, or any others)

8. The final averaged RDF data are in rdf.dat file, Plot the data with RDF vs. radial distance r.