Re: [AMBER] # of Coordinates Mismatch in PCA

From: Robert Molt <rwmolt07.gmail.com>
Date: Mon, 20 Aug 2018 16:50:18 -0400

Good afternoon,

It seems I have been two steps forward, one step back on this problem. I
am still encountering an error in matching the dimensionality of the
matrices. I believe that my masks are now perfectly consistent with one
another, unlike the preceding error. The error occurs only in the last
step of the below-listed ptraj script. The script is written in blue,
the ptraj output written in green.

parm piggy.prmtop
trajin full_no_waters.mdcrd 1 last 100
reference Crystal_Structure_with_Hydrogens.pdb
createcrd full_no_waters.mdcrd
run
crdaction full_no_waters.mdcrd rms ref
Crystal_Structure_with_Hydrogens.pdb :38-318&!.H=
crdaction full_no_waters.mdcrd matrix covar name full_no_waters_covar
:38-318&!.H=
runanalysis diagmatrix full_no_waters_covar out full_no_waters_evecs.dat
vecs 3 name myEvecs nmwiz nmwizvecs 3 nmwizfile protein.nmd nmwizmask
:38-318&!.H=
crdaction full_no_waters.mdcrd projection full_no_waters_projection3
modes myEvecs beg 1 end 1 :38-318&!.H=
clear all
readdata full_no_waters_evecs.dat name Evecs
parm piggy.prmtop
parmstrip !(:38-318&!.H=)

Up to this point, ptraj output clearly specifies that there are 1955
atoms after the mask, 2476 before. There are no errors. Here is the
specific text justifying my words:

   [parmstrip !(:38-318&!.H=)]
     Stripping atoms in mask [!(:38-318&!.H*)] (2476) from piggy.prmtop
     Stripped parm: 1955 atoms, 244 res, box: Orthogonal, 1 mol

parmwrite out cpu-gpu-mode1.prmtop
runanalysis modes name Evecs trajout cpu-gpu-mode1.nc \
   pcmin -100 pcmax 100 tmode 1 trajoutmask :38-318&!.H= trajoutfmt netcdf

At this point, I get the error:

   [runanalysis modes name Evecs trajout cpu-gpu-mode1.nc pcmin -100
pcmax 100 tmode 1 trajoutmask :38-318&!.H= trajoutfmt netcdf]
     Mask [:38-318&!.H*] corresponds to 1653 atoms.
     Writing 'cpu-gpu-mode1.nc' as Amber NetCDF
     ANALYZE MODES: Calculating coordinate projection using modes from Evecs
     Creating trajectory for mode 1
     Writing to trajectory cpu-gpu-mode1.nc
     PC range: -100.000000 to 100.000000
     Scaling factor: 1.000000
     Modes 'Evecs'
Error: # selected coords (4959) != eigenvector size (5865)

Now, the mask corresponds to 1653 atoms. I am confused; I have loaded no
new parameter files, and the mask is identical in each case. Seemingly,
though, the mask corresponds to a different number of atoms now. Can
someone point out the flaw in my logic?


On 8/14/18 7:55 PM, Jason Swails wrote:
>
>
>> On Aug 12, 2018, at 1:04 AM, Robert Molt <rwmolt07.gmail.com> wrote:
>>
>> Good morning,
>>
>> I am attempting to perform a principle component analysis but am
>> encountering an unexpected error during the projection step using cpptraj:
>>
>> Error: number selected coords (5865) != number avg coords (10983) in myEvecs
>>
>> In short, I am confused how a mismatch could occur. My understanding is
>> that if my eigenvectors (myEvecs) come directly from the covariance
>> matrix of a given trajectory, no mismatch should be possible (if I have
>> consistent masks). I must be confused and seek clarification.
>>
>> I am using cpptraj v. 15. I am following the tutorial specified on
>>
>> http://www.amber.utah.edu/AMBER-workshop/London-2015/pca/
>>
>> My cpptraj commands are (in blue), with trivial warning messages listed
>> (in red):
>>
>> trajin full_no_waters.mdcrd
>> reference Crystal_Structure_with_Hydrogens.pdb
>>
>> Warning: PDB line length is short (67 chars, expected 80).
>> Warning: PDB Crystal_Structure_with_Hydrogens.pdb: Reading frame 2, got
>> 711 atoms, expected 4431.
>> Warning: Only using frames 1-1
>>
>> createcrd full_no_waters.mdcrd
>> run
>> crdaction full_no_waters.mdcrd rms ref
>> Crystal_Structure_with_Hydrogens.pdb :38-318&!.H
> Note the mask. :38-318&!.H
>
> This selects all atoms in residues 38-318 not including atoms names H.
>
>> crdaction full_no_waters.mdcrd projection full_no_waters_projection3
>> modes myEvecs beg 1 end 3 :38-318&!.H=
> Note this mask with a = at the end. The = is a wildcard that omits all atoms that STARTS with H. This includes a lot more atoms that are excluded from the mask you’re doing the projection on.
>
>> Error: number selected coords (5865) != number avg coords (10983) in myEvecs
> Based on this, there are (10983-5865)/3=1706 atoms that start with H but are not named H.
>
> Make sure you’re using exactly the same mask in each step.
>
> HTH,
> Jason
>
> --
> Jason M. Swails
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber

-- 
Dr. Robert Molt Jr.
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Aug 20 2018 - 14:00:02 PDT
Custom Search