Why is the radius showing up in the fraction when the vectors are supposed to be normalized already? It's skewing up the direction of max increasing angle for the middle atom in the direction of the shortest side. What's the interpretation behind it?

PS: This also appears to be present in torsion calculations and oop angle gradients.

Hello
Thanks for the great youtube courses and utilities. I watched a many :)

I think the weights should be that of the most abundant isotope or relevant isotope in question, not the average weights. This is how other packages do it for ab-initio.

I use this reference:
https://physics.nist.gov/cgi-bin/Compositions/stand_alone.pl?ele=&ascii=ascii

I noticed while testing my C++ Z matrix code v psi4 which agree, but yours code gives very slightly different centre of mass coordinates when I add the following code (to your utility) to get the centre of mass adjusted Cartesian coordinates. Ultimately, this is caused by the weights you use.

Sorry. I don't know python to write code in normally, but this works okay.

print xyz coordinates to screen

def print_coords(self, comment):
    angstrom_to_bohr = 1.889725989
    total_mass = 0.0
    
    for i in range(self.n_atoms):
        total_mass += at_masses[self.atoms[i].attype]

    x_com = 0.0
    y_com = 0.0
    z_com = 0.0

    for i in range(self.n_atoms):
        x_com += self.atoms[i].coords[0] * at_masses[self.atoms[i].attype]
        y_com += self.atoms[i].coords[1] * at_masses[self.atoms[i].attype]
        z_com += self.atoms[i].coords[2] * at_masses[self.atoms[i].attype]
    
    x_com /= total_mass
    y_com /= total_mass
    z_com /= total_mass

    for i in range(self.n_atoms):
        self.atoms[i].coords[0] -= x_com
        self.atoms[i].coords[1] -= y_com
        self.atoms[i].coords[2] -= z_com
        
    print('\n')
    print('%s\n' % (comment), end='')
    for i in range(self.n_atoms):
        if self.atoms[i].attype == 'X':
                continue
        
        print('%-2s' % (self.atoms[i].attype), end='')
        
        for j in range(3):
            print(' %14.10f' % (self.atoms[i].coords[j] * angstrom_to_bohr), end='')
        
        print('\n', end='')

I think a bit higher precision in the printed output and definitions of weights would also be useful due to getting rounding errors you get from the cartesians you produce when used in calculations when comparing with others.

o' course weight are not used anyway in the utility to convert, so it is minor. The coordinates without COM adjustment are correct up to that point to within rounding errors.

Best regards
...

Hi tmpchem

I'm wondering what is the sorting convention you are using to align molecules with respect to the principal axis (ascending or descending)?. I'm trying to compute the gyration tensor of a macromolecule knowing beforehand its relationship to the inertia tensor (see appendix of dx.doi.org/10.1021/jp2065612).
In my case, all atoms have the same mass, and I would expect to get the same coordinates for both, the gyration tensor aligned molecule and inertia tensor aligned molecule.

Many thanks in advance for your time

I was looking at the calculations for ensemble properties of the system, and I found something confusing in regarding the calculation for kinetic pressure, defined as follows:

def get_virial(mol):
    """Clausius virial function for all atoms, force,s and coordinates.
    
    Args:
        mol (mmlib.molecule.Molecule): Molecule object with coordinate
            and force data.
    """
    mol.virial = 0.0
    for i in range(mol.n_atoms):
        for j in range(3):
            mol.virial += -mol.atoms[i].coords[j] * mol.g_total[i][j]

def get_pressure(mol):
    """Update total pressure of a system of molecules with boundary.
    
    Args:
        mol (mmlib.molecule.Molecule): Molecule object with temperature,
            volume, and virial data.
    """
    get_virial(mol)
    pv = mol.n_atoms * energy.kb() * mol.temp
    pv += mol.virial / (3.0 * mol.n_atoms)
    mol.press = kcalamol2pa() * pv / mol.vol

In equation form:

After spending an afternoon going through as much MD literature I could find, I've yet to see any paper that lists that N dividing the virial when calculating the pressure. Some examples:

Page 2: https://arxiv.org/pdf/1111.2705.pdf
Page 1: https://spiral.imperial.ac.uk/bitstream/10044/1/262/1/edm2006jcp.pdf
Page 5: http://phys.ubbcluj.ro/~tbeu/MD/C2_for.pdf

Not only that, but the literature is consistent as expected, and apparently in disagreement with this implementation. I am probably missing something, since I'm no expert in MD and the notation could be the culprit here. In any case, I'd like to know what's the problem.

when i run below line in jupyter notebook
run ./scripts/molecular_mechanics/mm.py ./geom/prm/co.prm

i get the following errors
ValueError Traceback (most recent call last)
~\computational_chemistry-master\scripts\molecular_mechanics\mm.py in
24
25 # Read in molecular geometry and topology.
---> 26 molecule = molecule.Molecule(input_file_name)
27
28 # Calculate potential energy.

~\computational_chemistry-master\scripts\molecular_mechanics\mmlib\molecule.py in init(self, infile_name)
480 self.GetTopology()
481 elif (self.filetype == 'prm'):
--> 482 self.ReadInPrm()
483 self.UpdateInternals()
484

~\computational_chemistry-master\scripts\molecular_mechanics\mmlib\molecule.py in ReadInPrm(self)
504 input_rows = fileio.GetFileStringArray(self.infile)
505
--> 506 self.atoms = fileio.GetAtomsFromPrm(input_rows)
507 self.bonds = fileio.GetBondsFromPrm(input_rows)
508 self.angles = fileio.GetAnglesFromPrm(input_rows)

~\computational_chemistry-master\scripts\molecular_mechanics\mmlib\fileio.py in GetAtomsFromPrm(rows)
401 for row in rows:
402 if row and row[0].upper() == 'ATOM':
--> 403 atoms.append(GetAtomFromPrm(row))
404 return atoms
405

~\computational_chemistry-master\scripts\molecular_mechanics\mmlib\fileio.py in GetAtomFromPrm(row)
231 if not _IsType(float, ro) or not float(ro) > 0.0:
232 raise ValueError(
--> 233 'Atomic vdw radius must be positive numeric value: %s' % ro)
234 if not _IsType(float, eps) or not float(eps) >= 0.0:
235 raise ValueError(

ValueError: Atomic vdw radius must be positive numeric value: 0.0