Author Topic: Partial Optimization (Intern Coord) of a big molecule  (Read 12979 times)

Arrepiadd

  • Jr. Member
  • **
  • Posts: 12
  • Karma: +0/-0
Partial Optimization (Intern Coord) of a big molecule
« on: April 23, 2008, 11:04:23 AM »
Hello everyone,

I am trying to optimize part of the structure of  a big system (94 atoms). I got the crystal structure so I got the positions of anything but H with a good accuracy. For the first steps I'd like to do an optimization of only the hydrogens.
I've been through the documentation and the tutorial and the only mention to fixing coordinates is in idef by using the command:
Code: [Select]
f stre x y
which fixed the bond length between atoms x and y. Then there are the bend and tors for the angles and dihedrals...

While this would be ok if I wanted to freeze a few atoms, I'd like to fix the positions of some 60 atoms. That adds up to a huge number of things to freeze and the possibility of doing something wrong should be quite high.

So, my question is is there a way to freeze everything for a list of atoms instead of going one by one? Something like:
Code: [Select]
f stre 1-60
that would freeze all the bonds that connect any atom between 1 and 60...
Because right now it seems that the easiest way to do it is editing the coord file by hand...

Thanks in advance.
Daniel

slackenerny

  • Jr. Member
  • **
  • Posts: 12
  • Karma: +0/-0
Re: Partial Optimization (Intern Coord) of a big molecule
« Reply #1 on: April 23, 2008, 06:15:39 PM »
IMHO you are doomed... :D
I solved a similar problem, by optimising carthesian instead of internal coordinates. That is quite demanding though it might be possible if you start with a reasonable startgeometry...
Basically either you invest a lot of your time or a lot of computational time. You decide what's worth more...


IMHO there should be a z2t tool (like x2t and t2x) to convert molden Z-matrices into a TM readable format... That would REALLY be a BIG held for problems like this and constrained geometries in general...

Arnim

  • Developers
  • Sr. Member
  • *
  • Posts: 238
  • Karma: +0/-0
Re: Partial Optimization (Intern Coord) of a big molecule
« Reply #2 on: April 23, 2008, 06:16:22 PM »
Hi!

If you want to optimize only the hydrogens, wouldn't it be easier to
simply fix the cartesian coordinates of the heavy atoms?
Just write f behind the atom
    0.00000000000000      0.00000000000000     -0.12178983933899      o f

Cheers,
Arnim

Arrepiadd

  • Jr. Member
  • **
  • Posts: 12
  • Karma: +0/-0
Re: Partial Optimization (Intern Coord) of a big molecule
« Reply #3 on: April 23, 2008, 06:56:35 PM »
Hi again!

Thanks to both for your comments.
Slackenerny's comment goes on the way of my line of thinking (especially the first and last lines  :) )

As for Arnim's comment, I thought of that already but that would lead to a cartesian optimization, right? Which, given the fact that I'm having SCF convergence problems with this molecule, will lead to an optimization taking... forever. So, if there is a way to put the f there and then asking for the internal redundants (which, I tried and could not put to work, but maybe it's my fault) then this is the thing I want. Otherwise, might as well just do the thing by hand because it would take me forever (hence, answering slackenerny's question: my time is less worth then computational time).

Once again, thanks for your help.

Best regards,
Daniel

uwe

  • Global Moderator
  • Sr. Member
  • *****
  • Posts: 487
  • Karma: +0/-0
Re: Partial Optimization (Intern Coord) of a big molecule
« Reply #4 on: April 23, 2008, 10:45:44 PM »
Hi,

as usual, there are more things possible with Turbomole than one might expect  ;)

1. yes, you can fix cartesian coordinates and optimize in redundant coordinates at the same time. This is an undocumented feature, and I am not sure if it is stable in your case and how efficient it will be, but if you fix the cartesian coords, define internal redundant coordinates (with ired in define), and then run the job with jobex -statpt [...], i.e. using the new statpt module instead of relax, it might work. At least with Turbomole 5.10

2. z-matrix for Turbomole? molden, babel, gaussview (via xyz format), or even the tmole script (which is part of Turbomole, see documentation chapter 1.9) can do that. But please remember that a z-matrix is a very intuitive method to define coordinates, but mathematically not really proper (one has to add dummy atoms to be able to describe a set of cartesian coordinates in some cases). Hence, using internal (redundant) coordinates is much more efficient than using z-matrices.

3. what is the reason for your SCF convergence problems? My experience with very difficult cases is that it might take a while to find the correct occupation and DIIS settings, etc. But once you have converged orbitals, a geometry optimization if often less demanding than the first energy calculation. Especially if only the hydrogens are allowed to move. Did you try $fermi ?

4. depending on which method you are using, the coordinates of experimental crystal structures do not have to be a good start structure for the optimization. The minimum geometry of HF, DFT, MP2, etc. can differ significantly from X-ray structures.

5. did you try to optimize the hole structure without constrains? RI-DFT with B-P86 functional, $marij and def-SV(P) basis set should not take too long for 94 atoms - but it depends if you have transition metal atoms, etc., so you might first try to get converged orbitals (see point 3).

Hope this helps a bit,

Uwe

Arrepiadd

  • Jr. Member
  • **
  • Posts: 12
  • Karma: +0/-0
Re: Partial Optimization (Intern Coord) of a big molecule
« Reply #5 on: April 24, 2008, 10:07:27 AM »
As always, Uwe comes and saves the day :)

as usual, there are more things possible with Turbomole than one might expect  ;)

Which is hard to find, given the scarce information on the documentation... :(

Quote
1. yes, you can fix cartesian coordinates and optimize in redundant coordinates at the same time. This is an undocumented feature, and I am not sure if it is stable in your case and how efficient it will be, but if you fix the cartesian coords, define internal redundant coordinates (with ired in define), and then run the job with jobex -statpt [...], i.e. using the new statpt module instead of relax, it might work. At least with Turbomole 5.10

Ok, we have turbomole 5.10 here so I tried this. I do not see any information anywhere on the fixed coordinates (not if I go into the 'i' menu in define or if I check the coordinates file). If I do it by hand, through 'idef' I see some fixed coordinates in a separate keyword ($intdef) in the coordinates file. So I may be doing something wrong. If this works as you say, then it's the thing I want. Didn't seem to work for me though...

Quote
2. z-matrix for Turbomole? molden, babel, gaussview (via xyz format), or even the tmole script (which is part of Turbomole, see documentation chapter 1.9) can do that. But please remember that a z-matrix is a very intuitive method to define coordinates, but mathematically not really proper (one has to add dummy atoms to be able to describe a set of cartesian coordinates in some cases). Hence, using internal (redundant) coordinates is much more efficient than using z-matrices.

Yes, I agree and I usually use internal redundant coordinates. But, if the description by z-matrices saves me a huge amount of time and allows me to do it more automatically, then in a few cases it may be worth it.

Quote
3. what is the reason for your SCF convergence problems? My experience with very difficult cases is that it might take a while to find the correct occupation and DIIS settings, etc. But once you have converged orbitals, a geometry optimization if often less demanding than the first energy calculation. Especially if only the hydrogens are allowed to move. Did you try $fermi ?

I haven't tried $fermi yet. I've been working with other software packages (which I am more familiar with) trying to get this calculation going but to no avail. So, I'm on my first few tries with Turbomole... the SCF problem comes from my experience with the other programs, not especially from the results of Turbomole usage. I'll give a look at $fermi though.

Quote
4. depending on which method you are using, the coordinates of experimental crystal structures do not have to be a good start structure for the optimization. The minimum geometry of HF, DFT, MP2, etc. can differ significantly from X-ray structures.

I know and ultimately the idea is to optimize the whole thing. But since we are going to compare the results with some solid stare NMR calculations and NMR is highly dependent on the geometry I want to use something as similar to what they have as possible. I'll also optimize the full geometry and see if it gets anything the experimentalists were not able to explain.

Quote
5. did you try to optimize the hole structure without constrains? RI-DFT with B-P86 functional, $marij and def-SV(P) basis set should not take too long for 94 atoms - but it depends if you have transition metal atoms, etc., so you might first try to get converged orbitals (see point 3).

I do not have transition metals (I have some Na+ though) and I did try that. I increased the SCF steps to 200 and the calculation stop after 200 SCF steps without converging...

As I said before, I went through the idef and fixed stuff by hand. I wanted to see the behavior of the program if I fixed stuff...
Since my idea is to fix anything that has (only) the atoms 1 through 60, it is just a matter of getting to the file and seeing which lines have only those numbers, so it shouldn't be hard to script.
All goes well if I look only and bond lengths. However, for angles and dihedrals the thing screws up a little bit... Here go the two sets of bonds, angles and dihedrals I fixed.
Code: [Select]
$intdef
# definitions of internal coordinates
   1 f  1.0000000000000 stre   18   25           val=   2.80071
   2 f  1.0000000000000 bend   14   16   13      val= 106.00700
   3 f  1.0000000000000 tors   10   14   13   12 val= -13.82794
   4 f  1.0000000000000 stre    6    5           val=   2.52608
   5 f  1.0000000000000 bend    3    1    2      val= 125.73701
   6 f  1.0000000000000 tors   45   44   41   40 val=   1.93098

And here is the same for the part automatically made by define (I'll change the numbering on the first column to the same numbers as above, just to make it more clear):
Code: [Select]
   1 k  1.0000000000000 stre   18   25           val=   2.80071
   2 k -0.6161601823403 bend   14   16   13      val=  20.30057
   3 k -0.8090169943750 tors   10   14   13   12 val=  25.24627
   3 k -0.5877852522925 tors   10   14   13   12 val= -18.09457
   4 k  1.0000000000000 stre    6    5           val=   2.52608
   5 k -0.5000000000000 bend    3    1    2      val=   3.16414
   5 k  0.8660254037844 bend    3    1    2      val=   4.79577
   5 k -1.0000000000000 bend    3    1    2      val=  -4.46234
   6 k -0.5000000000000 tors   45   44   41   40 val=  -8.27260
   6 k  0.8660254037844 tors   45   44   41   40 val= -12.92675
   6 k  1.0000000000000 tors   45   44   41   40 val=   8.54507

So, as you see, the 'val' is totally different! The manual says the internal coordinates do not have a direct meaning and cannot be changed by the user (page 46, first paragraph) but I do not want to change them, I want to fix them! Is there a way to get from the bottom information to the top one? Of course if I get your first suggestion to work, this is totally unnecessary...

Comments from anyone?

Thank you very much for your help, Uwe.

Best regards,
Daniel

uwe

  • Global Moderator
  • Sr. Member
  • *****
  • Posts: 487
  • Karma: +0/-0
Re: Partial Optimization (Intern Coord) of a big molecule
« Reply #6 on: April 25, 2008, 04:54:51 PM »
Hi,

Quote
Ok, we have turbomole 5.10 here so I tried this. I do not see any information anywhere on the fixed coordinates (not if I go into the 'i' menu in define or if I check the coordinates file). If I do it by hand, through 'idef' I see some fixed coordinates in a separate keyword ($intdef) in the coordinates file. So I may be doing something wrong. If this works as you say, then it's the thing I want. Didn't seem to work for me though...

Did you read the chapter about constrained optimizations (chapter 7.4) in the Tutorial?

First set the fixed internal coordinates with 'idef' and then 'f stre 1 2', etc. After that call ired to use internal redundant coordinates, including the frozen internal coordinates you have just defined.

If you look in the coord file, check the section

$redundant
     number_of_atoms            12
     degrees_of_freedom         30
     internal_coordinates       32
     frozen_coordinates          2
   Values of frozen coordinates
          2.62524980
          0.52359878
# definitions of redundant internals
   1 k  1.0000000000000 stre    1    7           val=   2.03717
   2 k -0.7071067811865 bend    2    7    1      val=   0.00000
        0.7071067811865 bend    6    7    1


As you can see, there are two frozen coordinates with the given values. The frozen internal coordinates are given at the end of the list:

  [...]
30 k -1.0000000000000 tors    2    3    4    5 val=   0.00000
        1.0000000000000 tors    3    4    5    6
       -1.0000000000000 tors    4    5    6    1
        1.0000000000000 tors    5    6    1    2
       -1.0000000000000 tors    6    1    2    3
        1.0000000000000 tors    1    2    3    4
  31 f  1.0000000000000 stre    1    2           val=   2.62525
  32 f  1.0000000000000 bend    3    2    1      val=  30.00000
[...]


Quote
Yes, I agree and I usually use internal redundant coordinates. But, if the description by z-matrices saves me a huge amount of time and allows me to do it more automatically, then in a few cases it may be worth it.

Just create your molecule with the z-matrix and save it as xyz. Convert the xyz file with 'x2t <filename> > coord' and say 'ired' in define. That's it.

Quote
  1 k  1.0000000000000 stre   18   25           val=   2.80071
   2 k -0.6161601823403 bend   14   16   13      val=  20.30057
   3 k -0.8090169943750 tors   10   14   13   12 val=  25.24627
   3 k -0.5877852522925 tors   10   14   13   12 val= -18.09457

O.k., this is a result from grep. I am afraid that this approach is too simple - the reason is that all those internal coordinates (except the first one) are linear combinations of primitive ones. Here a list of a complete set from a completely different molecule (so do not compare...):

Code: [Select]
  23 k  1.0000000000000 stre    5    6           val=   2.62525
  24 k  1.0000000000000 stre    6    1           val=   2.62525
  25 k -0.5000000000000 bend    1    3    2      val=   0.00000
       -0.5000000000000 bend    2    4    3
        1.0000000000000 bend    3    5    4
       -0.5000000000000 bend    4    6    5
       -0.5000000000000 bend    5    1    6
        1.0000000000000 bend    6    2    1

23 k and 24 k are simple bonds, but 25 k is a linear combination of 6 angles, and the coefficients of the linear combination are given in the third column (-0.5, -0.5, 1.0, ...). The val-value given here is that of the complete linear combination and has nothing to do with a 'real' angle.

This is a bit confusing for us human beings - I personally do not really know how such a linear combination of internal coordinates looks like or what the physical meaning of that combined internal coordinate is.  :)
Internal redundant coordinates are constructed the way that many, many (reasonable) internal coordinates are created (bond lengths, angles, torsions, out-of-plane angles, ...), and then the Wilson B-matrix is being diagonalized. The largest N eigenvalues and eigenvectors are taken as a complete set of internal redundant coordinates, because due to the larger eigenvalues, they are a mathematically more stable description of your coordinates than the z-matrix or the non-redundant internal coordinates. However, that makes them less intuitively understandable for us (or at least for me).

I would highly recommend to fix the cartesian coordinates of the non-hydrogens for your purpose, since you can only fix the complete linear combination of primitive internal coordinates, and not just one part of it.

Regards,

Uwe

Arrepiadd

  • Jr. Member
  • **
  • Posts: 12
  • Karma: +0/-0
Re: Partial Optimization (Intern Coord) of a big molecule
« Reply #7 on: April 25, 2008, 06:14:12 PM »
Hi,

Quote
Ok, we have turbomole 5.10 here so I tried this. I do not see any information anywhere on the fixed coordinates (not if I go into the 'i' menu in define or if I check the coordinates file). If I do it by hand, through 'idef' I see some fixed coordinates in a separate keyword ($intdef) in the coordinates file. So I may be doing something wrong. If this works as you say, then it's the thing I want. Didn't seem to work for me though...

Did you read the chapter about constrained optimizations (chapter 7.4) in the Tutorial?

Yes! What I meant with I tried and couldn't make it work is fixing the cartesian coordinates and then having the optimization in internal redundants.
I go into the coordinate file and add a 'f' to the end of each of the lines I want to fix (which are in cartesian coordinates).
I then go into define and ask for 'ired' and in the end of define I don't get the fixed coordinates (because inside the coordinates file I never get the frozen_coordinates up from 0).
Since, as you say, this is an undocumented feature, there is not much that I can do other than bother you :) Sorry for that.

Quote
First set the fixed internal coordinates with 'idef' and then 'f stre 1 2', etc. After that call ired to use internal redundant coordinates, including the frozen internal coordinates you have just defined.

Yes, that I know I can do and that is how I got the fixed coordinates I posted on my previous post. However, that process is the one I want to avoid. Fixing everything for 60 atoms is prone to errors. It works fine in that case, in which you fix one atom. To fix most of a (relatively) big structure it is not the best way. Too prone to errors...

Quote
Just create your molecule with the z-matrix and save it as xyz. Convert the xyz file with 'x2t <filename> > coord' and say 'ired' in define. That's it.

And can I fix distances in any way? By putting an f in the end of the line? By putting the values directly in the z-matrix instead of the end? By separating them into two parts?
No need to answer this, if it works, I prefer the internal redundants...

Quote
O.k., this is a result from grep. I am afraid that this approach is too simple - the reason is that all those internal coordinates (except the first one) are linear combinations of primitive ones. Here a list of a complete set from a completely different molecule (so do not compare...)
(...)
This is a bit confusing for us human beings - I personally do not really know how such a linear combination of internal coordinates looks like or what the physical meaning of that combined internal coordinate is.  :)
Internal redundant coordinates are constructed the way that many, many (reasonable) internal coordinates are created (bond lengths, angles, torsions, out-of-plane angles, ...), and then the Wilson B-matrix is being diagonalized. The largest N eigenvalues and eigenvectors are taken as a complete set of internal redundant coordinates, because due to the larger eigenvalues, they are a mathematically more stable description of your coordinates than the z-matrix or the non-redundant internal coordinates. However, that makes them less intuitively understandable for us (or at least for me).

Ok, so going into the file directly isn't an alternative. Thanks for the explanation.

Quote
I would highly recommend to fix the cartesian coordinates of the non-hydrogens for your purpose, since you can only fix the complete linear combination of primitive internal coordinates, and not just one part of it.

Yes, according to what you said two posts ago, that seems to be the thing I want. Fix the coordinates in cartesians and let the program make the fixed internal redundants based on that information. Unfortunately though, it is also the thing I cannot do...

Once again, thank you very much for your help. I am usually able to find the solution for my problems by myself (either by trial and error or by searching online) but with Turbomole I'm quite at loss...

Cheers,
Daniel

uwe

  • Global Moderator
  • Sr. Member
  • *****
  • Posts: 487
  • Karma: +0/-0
Re: Partial Optimization (Intern Coord) of a big molecule
« Reply #8 on: May 09, 2008, 05:22:57 PM »
Hi,

Quote
I go into the coordinate file and add a 'f' to the end of each of the lines I want to fix (which are in cartesian coordinates).
I then go into define and ask for 'ired' and in the end of define I don't get the fixed coordinates (because inside the coordinates file I never get the frozen_coordinates up from 0).

I see. The number of frozen coordinates in the $redundant section is only related to the internal redundant coordinates themselves. So if you fix cartesian coordinates, the value in the $redundant section will not be changed. Nevertheless, running a geometry optimization using internal redundant coordinates will keep the fixed cartesians at the same position.

To check that, just look at the cartesian coordinates in the file gradient - in gradient all coordinates and gradients of all steps of the optimization are saved, so you can compare the numbers.

CAUTION: You have to use statpt for the optimization, so call e.g.:

jobex -c 100 -ri -statpt

Without -statpt, jobex will use the old relax module which is not able to fix cartesian coordinates while optimizing in internal coodinates!

Another note: statpt (and therefore jobex) checks for the convergence of the internal redundant coordinates, and not of the cartesian ones. So while the geometry might be converged with respect to the gradient norm and energy, jobex will not recognize that.

You can change the convergence settings for statpt with $statpt (see chapter 12.2.15 of the Turbomole 5.10 documentation), e.g. by giving large thresholds for the displacements, but not for the gradient norm. See the job.last output and search for the section:

 
Quote
     ******************************************************************
                          CONVERGENCE INFORMATION

                               Converged?     Value      Criterion
             Energy change         no       0.0000106   0.0000010
             RMS of displacement   no       0.0368169   0.0005000
             RMS of gradient       yes      0.0003073   0.0005000
             MAX displacement      no       0.0982612   0.0010000
             MAX gradient          yes      0.0008484   0.0010000
      ******************************************************************

Let the job run (a lot of ?) cycles and check which ones of the convergence criteria are making problems and which ones are reasonable.

But I have no experience with such kinds of calculations, so perhaps it works quite well, perhaps it doesn't. It would be nice if you could give us some feedback after playing around with it.

Regards,

Uwe