Fixed Dihedral Optimisation help

[johal345]

Hi,

I’m using version 1.8 psithon and trying to fix all the dihedrals in my molecule so only the bond lengths can relax but I’m not sure what the correct syntax would be to do so.

Currently I’m using:

set frozen_dihedral 

but I’m not sure what to put to fix all the angles/ is there an options where I can just put * * * * to fix them all?

Thanks for any help!
Jay

[philipmnel]

Unfortunately there doesn’t seem to be any option to freeze every dihedral in a system. You’ll need to specify each to be frozen. You can pass frozen_dihedral a list of dihedrals like

set optking {
  frozen_dihedral = ("
    1 2 3 4
    5 6 7 8
  ")
}

[johal345]

Thanks for the reply, I came to the same conclusion looking through so have been testing the below to get my list:

import itertools
all_dihedrals = " ".join([str(atom) for dihedral in list(itertools.combinations(list(range(1,moleculename.natom()+1)), 4)) for atom in dihedral])

Is the new line important between the dihedral sets?

Cheers,
Jay

[philipmnel]

I don’t think the new line is important, just for readability.

[Rollin_King]

The downside of adding all possible number combinations is that the frozen coordinate will get added to the coordinate list. Since it’s frozen, this may not matter. If you wanted to avoid that, and if you don’t want (or can’t easily) add them by hand, you can get the list of all torsion coordinates that are generated by optking using the following:

params = p4util.prepare_options_for_modules()
optimizer_params = {k: v.get('value') for k, v in params.pop("OPTKING").items() if v.get('has_changed')}
opt_object = optking.opt_helper.CustomHelper(hab, params=optimizer_params)
intcos = opt_object.molsys.fragments[0].to_dict()['intcos']
for i in intcos:
if i['type'] == 'Tors':
     print(i['atoms'])

This is a quick and ugly hack, but it seems to work. Note it returns tuples of atom indices, that are zero-indexed. So you would need to add 1 to each number, compile them all in a string for your input, and it should work.

[johal345]

Thanks for this it makes a much more manageable list of torsions. Changed it slightly as don’t think I have the p4util package but seems to do the trick.

import optking
opt_object = optking.opt_helper.CustomHelper($moleculename$)
all_dihedrals = ""
for dct in opt_object.molsys.fragments[0].to_dict()['intcos']:
    if dct["type"] == "Tors":
        list_torsions = list(dct["atoms"])
        all_dihedrals += f"{list_torsions[0]+1} {list_torsions[1]+1} {list_torsions[2]+1} {list_torsions[3]+1} "