tmpchem / computational_chemistry

Files used in TMP Chem videos on computational chemistry
215 stars 93 forks source link

Z matrix to Cartesian converter uses average atomic weights #10

Open AlexB67 opened 4 years ago

AlexB67 commented 4 years ago

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 ...