[AMBER] PC modes per residue from existing evecs file?

From: Kenneth Huang <khuang8.student.gsu.edu>
Date: Mon, 7 May 2018 15:15:27 +0000

Hi all,


Is it possible to calculate the PC mode values on a per residue basis from an evecs file calculated from the entire system? Ie, if I calculate my original evecs file with


parm example.prmtop

trajin example.nc 1 last 1

rms first :4-101&!H=

matrix name matc covar :4-101&!H=
run
runanalysis diagmatrix matc out PCA.evecs.dat name diagmodes vecs 5
rms first :4-101&!H=
projection evecs diagmodes out PCA.project.dat beg 1 end 5 :4-101&!H=

And then in essence, want to construct a projection but only for a single residue on the same trajectory in something like-

parm example.prmtop

trajin example.nc 1 last 1

readdata test.dat name diagmodes

rms first :4-101&!H=

projection evecs diagmodes out new.dat beg 1 end 5 :4&!H=

I know the above is incorrect, and won't work because projection wants the same atom selection in matrix since the header on the evecs file is different, etc. I've had mild success in manipulating the evecs file by removing pieces until I have the matching entries corresponding to that residue, but the evecs file formatting isn't the easiest to script around.


So I was wondering, is there an approximate way to backtrack from a calculated mode, say PC1, and deconstruct it back into its component parts on a per residue basis? I've dug around the source code and I'm guessing the relevant part is


    for (int mode = beg_; mode < end_; ++mode) {

      DataSet_Modes::AvgIt Avg = modinfo_->AvgBegin();

      double proj = 0;

      std::vector<double>::const_iterator sqrtmass = sqrtmasses_.begin();

      for (AtomMask::const_iterator atom = mask_.begin(); atom != mask_.end(); ++atom)

      {

        const double* XYZ = frm.Frm().XYZ( *atom );

        double mass = *(sqrtmass++);

        proj += (XYZ[0] - *(Avg++)) * mass * Vec[0];

        proj += (XYZ[1] - *(Avg++)) * mass * Vec[1];

        proj += (XYZ[2] - *(Avg++)) * mass * Vec[2];

        Vec += 3;

      }

      float fproj = (float)proj;

      project_[mode]->Add( frameNum, &fproj );

    }

Which I'm reading very crudely as eigenvector (x,y,z components) times atom mass times average position (xyz summed?) times the eigenvalue constant?


Best,

Kenneth
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon May 07 2018 - 08:30:02 PDT
Custom Search