getting the center of mass of each part of a molecule

[jladasky at itu.edu]

On Saturday, May 20, 2017 at 9:13:39 AM UTC-7, qasi... at gmail.com wrote
> I have more than 100 ligand molecules. The image I showed at the beginning of this discussion was only an example. It had missing atoms. 
> See it without missing any atom in the pdb format (coordinates in Angstrom unit):

[snip atomic coordinate tables]

> I am doing molecular dynamics simulation of protein-ligand complex. The issue is to keep the ligand in the binding site of the protein during simulation when turned off the interactions (van der waals, coulomb...) between the ligand and protein. In that case, to keep the ligand in the bindind site I need to apply the translational and rotational restraints (distance, angle, dihedral restraints). The more detail is here http://pubs.acs.org/doi/abs/10.1021/jp0217839

Well now, that's Ph.D level work.  I work with GROMACS myself.  

I tried to obtain a copy of the paper.  Why is a 14 year old article behind a paywall?  I requested a reprint from the authors through ResearchGate.
 
> My aim is to get the COM of the whole ligand (I already do it) and divide the ligand in two parts and find to the COM of each part. Then I will apply the distance, angle and dihedral restraints on the atoms of the ligand and protein.  
> -The first distance restraints will be between one of the ligand heavy atom closest to the ligand COM (COM for whole ligand atoms) and one of the protein atoms closest to the ligand COM. 
> -The second distance restraints  will be between one of the ligand heavy atom closest to the COM of the first part of ligand in two parts and one of the protein atoms closest to the COM of the first part of ligand in two parts. 
> -And the third distance restraints  will be between one of the ligand heavy atom closest to the COM of the second part of ligand in two parts and one of the protein atoms closest to the COM of the second part of ligand in two parts.
> 
> That is, I will form 3 bonds between ligand and protein. Hope that makes clear what I am trying to do.

Is that the method that was used in the research paper, or is that something that you devised on your own?
 
> Please tell me on the python code what I need to do.

Your requirements are still rather more vague than you realize.  I think that you are counting on the fact that the atomic coordinates in a PDB file include a certain amount of noise, and possibly also any strain in the ligand molecule.  If you have a truly symmetrical ligand, your center of mass computations will not be robust.

But I see the general idea behind what you are trying to do.  And I suppose any method that forms three virtual bonds between three distinct atoms in ligand and the protein will do the job.  Maybe you don't care that you get the true center of mass, just that you get close to it.

How much Python do you want to learn?  I am going to suggest a basic modification to your approach to coding.  You will have several steps in this algorithm.  Learn to structure your code in in a way which allows you to read your algorithm in plain language using comments, descriptive variable names, or descriptive function names.  

Right now you show one function called heavy_atoms.  That is not a very descriptive name.  I'm not even sure why it's a function at all, because it doesn't return any values to whatever called it.  If your program is short and you only want to work with global variables, you can skip the use of functions.

I would like to see your whole program.  However, at this stage it's probably not of interest to other readers of the newsgroup.  Please consider continuing through email.  Thanks.