The current version of the script assumes that there is only one glycerol backbone in the lipid. This is, however, not the case e.g. for cardiolipin.
The assumption does not really matter if there are any lipids with just one glycerol backbone (PC, PG, PE, ...) in the bilayer. However, if the bilayer would be of pure cardiolipin, the centering script would not work. It would (unfortunately) not crash, but instead of at step 2/3 centering around one of the tail methyl groups, it would center around whatever atom happens to be last in the concatenated mapping files. Typically this will be water or an ion, which will mess up the whole centering.
The fix: Script that finds also the glycerols of lipids that have multiple glycerols.
localhost:POPCandCLinCHARMM36 bb$ diff -y example.orig example.new
Mapping: mappingMolecules.txt Mapping: mappingMolecules.txt
TPR in: PC_CL.tpr TPR in: PC_CL.tpr
XTC in: short.xtc XTC in: short.xtc
XTC out: out.xtc XTC out: out.xtc
Will now CENTER trajectory short.xtc to out.xtc. Will now CENTER trajectory short.xtc to out.xtc.
Making molecules whole... Making molecules whole...
Selected 0: 'System' Selected 0: 'System'
Reading frame 0 time 0.000 Reading frame 0 time 0.000
Reading frame 20 time 400.000 -> frame 20 time Reading frame 20 time 400.000 -> frame 20 time
The name of the CH3 carbon in sn-1 chain: C316 | The name of the CH3 carbon in sn-1 chain: CD18
The atom number of the last lipid's CH3: 22917 | Found 60 atoms with name CD18
Centering molecules around atom 22917... | The atom number of the last lipid's CH3: 14877
> Centering molecules around atom 14877...
Selected 7: 'centralAtom' Selected 7: 'centralAtom'
Selected 0: 'System' Selected 0: 'System'
The name of the g3 carbon(s) in POPC: C1 The name of the g3 carbon(s) in POPC: C1
No g3 carbons in TIP3. No g3 carbons in TIP3.
The name of the g3 carbon(s) in TOCL2: C31 C11 The name of the g3 carbon(s) in TOCL2: C31 C11
No g3 carbons in SOD. No g3 carbons in SOD.
Centering around the center of mass of index group g3carbons. Centering around the center of mass of index group g3carbons.
Selected 8: 'g3carbons' Selected 8: 'g3carbons'
Selected 0: 'System' Selected 0: 'System'
Centering done Centering done
Calculating density profiles... Calculating density profiles...
rm: electronDENSITY.xvg: No such file or directory rm: electronDENSITY.xvg: No such file or directory
- POPC - - POPC -
The groups selected for centering and for density calculation The groups selected for centering and for density calculation
Selected 8: 'g3carbons' Selected 8: 'g3carbons'
Selected 3: 'POPC' Selected 3: 'POPC'
- TIP3 - - TIP3 -
The groups selected for centering and for density calculation The groups selected for centering and for density calculation
Selected 8: 'g3carbons' Selected 8: 'g3carbons'
Selected 4: 'TIP3' Selected 4: 'TIP3'
- TOCL2 - - TOCL2 -
The groups selected for centering and for density calculation The groups selected for centering and for density calculation
Selected 8: 'g3carbons' Selected 8: 'g3carbons'
Selected 2: 'TOCL2' Selected 2: 'TOCL2'
- SOD - - SOD -
The groups selected for centering and for density calculation The groups selected for centering and for density calculation
Selected 8: 'g3carbons' Selected 8: 'g3carbons'
Selected 5: 'SOD' Selected 5: 'SOD'
Calculating form factor... Calculating form factor...
which shows that the fixed code picks a cardiolipin tail CH3-carbon instead of a POPC one (because TOCL2 was later in the mapping file than POPC). Despite of this, however, the outcome is practically unchanged:
The current version of the script assumes that there is only one glycerol backbone in the lipid. This is, however, not the case e.g. for cardiolipin.
The assumption does not really matter if there are any lipids with just one glycerol backbone (PC, PG, PE, ...) in the bilayer. However, if the bilayer would be of pure cardiolipin, the centering script would not work. It would (unfortunately) not crash, but instead of at step 2/3 centering around one of the tail methyl groups, it would center around whatever atom happens to be last in the concatenated mapping files. Typically this will be water or an ion, which will mess up the whole centering.
The fix: Script that finds also the glycerols of lipids that have multiple glycerols.
Demonstration
Run original:
Fix the code. Run the fixed code:
Compare:
which shows that the fixed code picks a cardiolipin tail CH3-carbon instead of a POPC one (because TOCL2 was later in the mapping file than POPC). Despite of this, however, the outcome is practically unchanged: