catdcd output and .psf file atom count mismatch

From: Izhar Karbat (izhar.karbat_at_gmail.com)
Date: Mon May 06 2013 - 17:18:53 CDT

Hi all.

I'm dealing with a long trajectory of a membrane embedded protein complex, divided between 3 separate DCD files. for the sake of analysis, i would like to extract only the protein coordinates( discarding the lipids/water/ions etc .) into a single DCD that would contain the entire trajectory. Reading previous posts on this subject i have used the following piece of code, which generates an index file containing the protein atom indexes (findexfile.ind), and a PSF containing these atoms only (complex_prot.psf ).

******************************
package require psfgen
topology ../../../top_all27_prot_lipid.rtf

  
readpsf ../complex_popcwi.psf
coordpdb ../complex_popcwi.pdb

mol load psf ../complex_popcwi.psf pdb ../complex_popcwi.pdb
  
set all [atomselect top all]
$all set beta 0
set sel1 [atomselect top "protein"]
  
set indices [$sel1 get index]
  
set file [open findexfile.ind w]
foreach i $indices {
            puts $file $i
            set sel [atomselect top "index $i"]
            $sel set beta 1
}
  
flush $file
close $file

set NotProtein [atomselect top "beta = 0"]

set seglist [$NotProtein get segid]
set reslist [$NotProtein get resid]

foreach segid $seglist resid $reslist
      delatom $segid $resid
}
  
writepsf complex_prot.psf
*******************************

The call to catdcd is as follows :

./catdcd -o protein_trajectory.dcd -i findexfile.ind complex_popcwieq-04.dcd complex_popcwieq-05.dcd complex_popcwieq-06.dcd

 the script prepares the two files at the very same loop, first writing the atom index to findexfile.ind, and then mark it for preservation for the PSF file.
Yet, when trying to load the resulting trajectory with the complex_prot.psf protein structure file, i get the following error :

ERROR) BaseMolecule: attempt to init atoms while structure building in progress!
ERROR) Invalid number of atoms in file: 7689

indeed, the .psf contains only 7677 atoms.

it is unclear to me what could be the source to this discrepancy (so basically i am stuck :( ).
any help would be highly appreciated.

This archive was generated by hypermail 2.1.6 : Wed Dec 31 2014 - 23:21:11 CST