The script failed in the case where lipids are made of residues, as in the Amber force field.
Here I fix this by asking the user to provide an index file (missingMolecules.ndx) if the molecules found in the molecule mapping file are not found in the tpr-file.
Another issue that emerged was that some groups (such as NA in the below example) were found more than once when creating the index from the tpr-file. I fixed this such that now the group numbers (instead of group names) are used in the gmx density step.
Demonstration
Failure of the original script:
localhost:CLinAmber bb$ pwd
MATCH/scratch/scriptsbyhsantila/testSystems/CLinAmber
localhost:CLinAmber bb$ ../../center_n_calcFF.sh mappingMolecules.txt CL.tpr short.xtc out.xtc
Mapping: mappingMolecules.txt
TPR in: CL.tpr
XTC in: short.xtc
XTC out: out.xtc
Will now CENTER trajectory short.xtc to out.xtc.
Making molecules whole...
Selected 0: 'System'
Reading frame 0 time 60000.000
Last frame 5 time 61000.000 -> frame 5 time 61000.000
The name of the CH3 carbon in sn-1 chain: CBI
Found 90 atoms with name CBI
The atom number of the last lipid's CH3: 22317
Centering molecules around atom 22317...
Selected 14: 'centralAtom'
Selected 0: 'System'
The name of the g3 carbon(s) in TOCL2: CAB2 CCD2
No g3 carbons in SOL.
No g3 carbons in NA.
Centering around the center of mass of index group g3carbons...
Selected 15: 'g3carbons'
Selected 0: 'System'
Centering done
Calculating density profiles...
rm: electronDENSITY.xvg: No such file or directory
- TOCL2 -
gmx density failed, exiting.
localhost:CLinAmber bb$
After fixing this such that missingMolecules.ndx is required, if the molecule is not found in the standard index file (created from the tpr-file), we get another fail:
localhost:CLinAmber bb$ ../../center_n_calcFF.sh mappingMolecules.txt CL.tpr short.xtc out.xtc
Mapping: mappingMolecules.txt
TPR in: CL.tpr
XTC in: short.xtc
XTC out: out.xtc
Will now CENTER trajectory short.xtc to out.xtc.
Making molecules whole...
Selected 0: 'System'
Reading frame 0 time 60000.000
Last frame 5 time 61000.000 -> frame 5 time 61000.000
The name of the CH3 carbon in sn-1 chain: CBI
Found 90 atoms with name CBI
The atom number of the last lipid's CH3: 22317
Centering molecules around atom 22317...
Selected 14: 'centralAtom'
Selected 0: 'System'
The name of the g3 carbon(s) in TOCL2: CAB2 CCD2
No g3 carbons in SOL.
No g3 carbons in NA.
Centering around the center of mass of index group g3carbons...
Selected 16: 'g3carbons'
Selected 0: 'System'
Centering done
Calculating density profiles...
- TOCL2 -
The groups selected for centering and for density calculation:
Selected 16: 'g3carbons'
Selected 15: 'TOCL2'
- SOL -
The groups selected for centering and for density calculation:
Selected 16: 'g3carbons'
Selected 6: 'SOL'
- NA -
gmx density failed:
[...]
Fatal error:
Cannot read from input
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------
Selected 16: 'g3carbons'
Error: Multiple groups 'NA' selected
Error: No such group 'NA'
localhost:CLinAmber bb$
To avoid the problem of several identical groups with identical names, let us use group numbers instead. (These are extracted from the foo.ndx by grepping.) Now the script runs to the end:
localhost:CLinAmber bb$ ../../center_n_calcFF.sh mappingMolecules.txt CL.tpr short.xtc out.xtc
Mapping: mappingMolecules.txt
TPR in: CL.tpr
XTC in: short.xtc
XTC out: out.xtc
Will now CENTER trajectory short.xtc to out.xtc.
Making molecules whole...
Selected 0: 'System'
Reading frame 0 time 60000.000
Last frame 5 time 61000.000 -> frame 5 time 61000.000
The name of the CH3 carbon in sn-1 chain: CBI
Found 90 atoms with name CBI
The atom number of the last lipid's CH3: 22317
Centering molecules around atom 22317...
Selected 14: 'centralAtom'
Selected 0: 'System'
The name of the g3 carbon(s) in TOCL2: CAB2 CCD2
No g3 carbons in SOL.
No g3 carbons in NA.
Centering around the center of mass of index group g3carbons...
Selected 16: 'g3carbons'
Selected 0: 'System'
Centering done
Calculating density profiles...
- TOCL2 -
The groups selected for centering and for density calculation:
Selected 16: 'g3carbons'
Selected 15: 'TOCL2'
- SOL -
The groups selected for centering and for density calculation:
Selected 16: 'g3carbons'
Selected 6: 'SOL'
- NA -
The groups selected for centering and for density calculation:
Selected 16: 'g3carbons'
Selected 4: 'NA'
Calculating form factor...
localhost:CLinAmber bb$
Also, it gives an output identical to the original version:
localhost:CLinAmber bb$ cd ../POPCandCLinCHARMM36/
localhost:POPCandCLinCHARMM36 bb$ ../../center_n_calcFF.sh mappingMolecules.txt PC_CL.tpr short.xtc out.xtc
Mapping: mappingMolecules.txt
TPR in: PC_CL.tpr
XTC in: short.xtc
XTC out: out.xtc
Will now CENTER trajectory short.xtc to out.xtc.
[...]
- SOD -
The groups selected for centering and for density calculation:
Selected 8: 'g3carbons'
Selected 5: 'SOD'
Calculating form factor...
localhost:POPCandCLinCHARMM36 bb$ diff Electron_Density_From_Simulation.dat orig/Electron_Density_From_Simulation.dat
localhost:POPCandCLinCHARMM36 bb$
The script failed in the case where lipids are made of residues, as in the Amber force field.
Here I fix this by asking the user to provide an index file (missingMolecules.ndx) if the molecules found in the molecule mapping file are not found in the tpr-file.
Another issue that emerged was that some groups (such as NA in the below example) were found more than once when creating the index from the tpr-file. I fixed this such that now the group numbers (instead of group names) are used in the gmx density step.
Demonstration
Failure of the original script:
After fixing this such that missingMolecules.ndx is required, if the molecule is not found in the standard index file (created from the tpr-file), we get another fail:
To avoid the problem of several identical groups with identical names, let us use group numbers instead. (These are extracted from the foo.ndx by grepping.) Now the script runs to the end:
Also, it gives an output identical to the original version: