Optimization with additional harmonic constraints

[bgobbi]

Hi Everybody,
I would like to run a geometry optimization with additional harmonic (quadratic) constraints.
The goal is to optimize the most strained part of an input conformation only. Aiming see which part of the molecule relaxes most and in which direction. Other QM packages have the same feature (Gaussian, Gamess).

Internally the geometry optimization adds a quadratic potential to the deviation from the input geometry.
eg.

- if atom1 is in its original x,y,z coordinate it gets a penalty of 0
- if atom2 is 0.1A from its original xyz position a penalty of
  k * (0.1)^2 is added to the total energy

k would be a user defined parameter.

Ideally this could be done on a whole molecule (all atom have the same k) or per atom basis (each atom could have a different k).

Does anybody know if this is doable in Psi4 or how difficult it would be to implement?

Thanks,
Alberto

[Rollin_King]

There is a mechanism in PSI4 to specify a 'target' value of an internal coordinate with fixed_distance(), fixed_bend(), and fixed_dihedral() in the input. These add harmonic force constants progressively larger around the specified values (the original intent was to optimize toward a specified value not already satisfied, but it could be the same). We need a workaround for cartesians and to make the value constant.

Here is something that might work:

1. Run an optimization with

   set optking {
   opt_coordinates = cartesian
   intcos_generate_exit true
   } # in order to generate an .intco file.

2. Modify the .intco file by adding the equilibrium value (around which the force constant will be applied) to the end of the desired line (e.g., to the end of "X 1 Z" for the z-coordinate of atom 1. This will be the same as the initial cartesian coordinate in your case.

3. In order to make the applied force constant a constant you will have to modify the source code slightly (see in psi4/optking/molecule MOLECULE::apply_constraint_forces()).

[yuhong11]

Hi Mr. King,

Could you please advise if this is the correct way of modifying the source as you mentioned?
I found the function void MOLECULE::apply_constraint_forces(void) in optking/molecule.cc file. In line 196, it shows,

k = (1 + 0.05 * (iter-1)) * Opt_params.fixed_coord_force_constant;

And will this work if I change it to

k = Opt_params.fixed_coord_force_constant;

I will compile psi4 after this change.
Thank you very much for your time.
Yuhong

I will paste the whole function down here.

// Apply extra forces for internal coordinates with user-defined
// equilibrium values.
void MOLECULE::apply_constraint_forces(void) {
double * f_q = p_Opt_data->g_forces_pointer();
double **H = p_Opt_data->g_H_pointer();
int N = Ncoord();
int iter = p_Opt_data->g_iteration();
double k;
int cnt = -1;
for (std::size_t f=0; f<fragments.size(); ++f) {
for (int i=0; i<fragments[f]->Ncoord(); ++i) {
++cnt;
if (fragments[f]->coord_has_fixed_eq_val(i)) {
double eq_val = fragments[f]->coord_fixed_eq_val(i);
double val = fragments[f]->coord_value(i);
// Increase force constant by 5% of initial value per iteration
k = (1 + 0.05 * (iter-1)) * Opt_params.fixed_coord_force_constant;
H[cnt][cnt] = k;
double force = (eq_val - val) * k;
oprintf_out("\tAdding user-defined constraint: Fragment %d; Coordinate %d:\n", f+1, i+1);
oprintf_out("\t\tValue = %12.6f; Fixed value = %12.6f\n", val, eq_val);
oprintf_out("\t\tForce = %12.6f; Force constant = %12.6f\n", force, k);
f_q[cnt] = force;
// If user eq. value is specified delete coupling between this coordinate and others.
oprintf_out("\tRemoving off-diagonal coupling between coordinate %d and others.\n", cnt+1);
for (int j=0; j<N; ++j)
if (j != cnt)
H[j][cnt] = H[cnt][j] = 0;
}
}
}
return;
}

[Rollin_King]

Yes, your change will allow the user to specify a force constant that will be used to add a force to any coordinate for which an 'equilibrium value' (in your case, the initial value) is given. We don't have a FIXED_CARTESIANS at present (which would be nice to add), which is why I suggested you try modifying the intco file by hand to add a force according to a cartesian displacement. I've not tried this, so I can't guarantee it will work.

The provided force constant is used by the approximate Hessian in the optimization (H[cnt][cnt] = k) which may or may not be advantageous to your optimization. If k is large, then this is desirable.

You will also need to change another line. 'f_q[cnt] = force' to 'f_q[cnt] += force' since you want to add your extra force to the existing ones.

[yuhong11]

Dear Dr. KIng,

I dig thru the source and found that
Opt_params.fixed_coord_force_constant ( which is k ) is assigned to 0.5 in these lines:
/psi4/optking/set_params.cc:423: Opt_params.fixed_coord_force_constant = options.get_double("FIXED_COORD_FORCE_CONSTANT");
./read_options.cc:2449: options.add_double("FIXED_COORD_FORCE_CONSTANT", 0.5);

But say if I need this k to be exactly 1.0, do I just set k = 1.0 and leave the rest unchanged.

Also do I need to prepare two input files to run optimization with cartesian constraints? B/c I need to first modify the .intco file. Then I rerun the input without set intcos_generate_exit true? Could you please provide me some input file examples?

Thank you very much,
Yuhong

[Rollin_King]

The line in read_options.cc is simply how we set our defaults for options. If a value is specified in the input, then this default value is not used. So you shouldn't need to change the defaults in the source code for any of the options.

Yes, since this feature is not built-in at present. You will need to either modify the input file, or construct two different ones. I have never attempted what you are doing with cartesians, but there is a reasonable chance of success.

[yuhong11]

Hi Dr. King,

Sorry for bothering you so much. Could you please explain what f_q[cnt] = force does?
Also could you tell me that how to set up the k value in the input if that's possible.

Thank you very much
Yuhong