Control and Prediction of the Organic Solid State

Using DMACRYS and NEIGHCRYS

The new programmes, NEIGHCRYS and DMACRYS, can be used for larger systems than were previsouly possible with DMAREL and Neighbours. This is because of the new numbering system used. The codes can also use a wider range of potentials - either user-defined or, more significantly, the Williams 01 potentials [Williams, D.E., J. Comp. Chem (2001)], where the implicite foreshortening of distances to hydrogen is dealt with.

Here is an outline of how to use the new codes in the following sections:

General notes on the new codes
Changes to the labelling and file contents (fort.21, fort.12, mol.geom, dmarel.axis, fit.pots, bondlengths)
How to do a proper Williams run (for a single res file)
How to do a ‘wrong’ Williams run
How to do an old FIT run

1. General notes on the new codes

The main feature of the new codes is that a greter number of atoms can be used - up to 99999 atoms in total may be present in the unit cell. This has necessitated a change in the atom labelling used, which means that all files are a different format.

The new utility programme NEIGHCRYS will correctly generate the .dmain file from your starting files.

DMACRYS is also able to calculate the properties using the same programme (for DMAREL, you used to use DMAREL3.11 for the properties calculation), and pressure can be introduced into the system.

2. Changes to the labelling and file contents (fort.21, fort.12, mol.geom, dmarel.axis, fit.pots, bondlengths)

Labels

The atom labelling system has changed since the last version of DMAREL was released.
The atoms are labelled C_F1_1___1___. The underscores are essential, and the label is broken down like this:
C_ : The first two characters in the label are the atom type. If there is only letter in the label, the blank is filled with an underscore. E.g. C_ ,O_ ,N_, H_, Br, Cl, etc. Note that these labels are case sensitive. This means that is the SFAC line of your .res file (where NEIGHCRYS reads the atoms types in) the first letter must be capital and the second small.
F1 : The next two characters are the potential type. These are assigned automatically, but you should check them manually. If you are using FIT, the types are all F1, except the hydrogens which may be F2 or F3.
_ : The next character is a flag. This is blank, unless the atom is in a moleucle that has been inverted by NEIGHCRYS, in which case it is I.
1____ : The next 5 characters are the atom number in the molecule. Numbering restarts for each atom type, i.e. all carbon atoms are labelled 1-5 (say), then nitrogen 1-2, oxygen 1-2 and hydrogen 1-8 (see pasted text below) This completes the short atom label.
1____ : The last 5 characters complete the full atom label. This identifies the atom number in the unit cell, and sets the limit of the maximum number of atoms in the cell to 99999. (see pasted text below)

These atom labels will infiltrate many of the files. You can find these labels in fort.21. To generate them, run NEIGHCRYS with only the basic files – your .res file, a potential file, and bondlengths. NB: The bondlengths file must have the atom labels changed from CAR to C_ etc.
This will generate a list of atom types in your fort.21 file that looks like this:

          Equivalent basis atoms
      atom atomic  name  input molecule invert
     index number        name   number   flag

          1 6 C_F1_1____1____  C1 1 F
          2 6 C_F1_1____15___     2 F
          3 6 C_F1_2____2____  C2 1 F
          4 6 C_F1_2____16___     2 F
          5 8 O_F1_1____3____  O1 1 F
          6 8 O_F1_1____17___     2 F
          7 8 O_F1_2____4____  O2 1 F
          8 8 O_F1_2____18___     2 F
          9 7 N_F1_1____5____  N1 1 F
          10 7 N_F1_1____19___    2 F
          11 7 N_F1_2____6____ N2 1 F
          12 7 N_F1_2____20___    2 F
          13 7 N_F1_3____7____ N3 1 F
          14 7 N_F1_3____21___    2 F
          15 7 N_F1_4____8____ N4 1 F
          16 7 N_F1_4____22___    2 F
          17 1 H_F2_1____9____ H1 1 F
          18 1 H_F2_1____23___    2 F
          19 1 H_F2_2____10___ H2 1 F
          20 1 H_F2_2____24___    2 F
          21 1 H_F2_3____11___ H3 1 F
          22 1 H_F2_3____25___    2 F
          23 1 H_F2_4____12___ H4 1 F
          24 1 H_F2_4____26___    2 F
          25 1 H_F2_5____13___ H5 1 F
          26 1 H_F2_5____27___    2 F
          27 1 H_F2_6____14___ H6 1 F
          28 1 H_F2_6____28___    2 F

Helpfully, this acts as a lookup table, with the atom label from the res file and the DMACRYS label.

dmarel.axis

Using DMAREL and Neighbours, this file looked like this:

        MOLX 1 
        X LINE CAR1 CODA CAR2 CODA 2
        Y PLANE CAR1 CODA CAR2 CODA 2 OXY1 CODA 2
        ENDS

But to use with DMACRYS and NEIGHCRYS, you should adjust it to use the new atom labels, and drop the CODA keyword:

MOLX 1                                                                          
X LINE   C_F1_2____ C_F1_3____ 1
Y PLANE  C_F1_2____ C_F1_3____ 1 C_F1_1____ 1                                   
ENDS

Notice that is uses the short atom labels. Also, the potential type must be correct.

Fort.21

With this axis file, you can now generate your molecular conformation for Gaussian. Run NEIGHCRYS, with the files as before and your new dmarel.axis file. Now, in fort.21 you will find the following familiar section:

Local axes from input dataset

Local axis set for molecule: 1
x=> -0.891615 -0.347799 -0.289928
y=> 0.199207 -0.876318 0.438615
z=> -0.406619 0.333320 0.850622

Atom positions in local axis system for molecule 1

basis No. Species x y z (Angstroms) Mass

1 C_F1_1____ -1.457148 -1.138491 -0.000054 12.010700
9 C_F1_2____ -0.776877 0.017477 -0.000054 12.010700
17 C_F1_3____ 0.695207 0.017477 -0.000054 12.010700
25 O_F1_1____ 1.380979 1.034138 0.000194 15.999400
33 H_F1_1____ -0.936056 -2.087754 -0.001291 1.007940
41 H_F1_2____ -2.536577 -1.170809 -0.000373 1.007940
49 H_F1_3____ -1.270829 0.980410 0.001203 1.007940
57 H_F1_4____ 1.159371 -0.987246 -0.000686 1.007940

basis No. Species x y z (AU) Mass

1 C_F1_1____ -2.753612 -2.151436 -0.000102 12.010700
9 C_F1_2____ -1.468086 0.033026 -0.000102 12.010700
17 C_F1_3____ 1.313752 0.033026 -0.000102 12.010700
25 O_F1_1____ 2.609672 1.954238 0.000367 15.999400
33 H_F1_1____ -1.768890 -3.945285 -0.002440 1.007940
41 H_F1_2____ -4.793437 -2.212509 -0.000705 1.007940
49 H_F1_3____ -2.401519 1.852707 0.002274 1.007940
57 H_F1_4____ 2.190895 -1.865626 -0.001297 1.007940

This can then be used by GAUSSIAN to generate the wavefunction, for use by NEIGHCRYS.

3. The new Williams Potential

We have all heard of the "Williams" potential as an update to our present "FIT" potential. Presently our atom-atom potentials define types based on atomic number (i.e CAR, NIT, OXY), with only one set of parameters per type. The exception to this is hydrogen, which has different parameters if it is bonded to oxygen, or nitrogen.
The Williams potential expands on this, with parameters for different bonding scenarios, i.e. 2-, 3- and 4-coordinated carbon atoms; NH₂⁺, -NH, NC₃ and N=N nitrogen, etc. This makes a lot of chemical sense, and should improve our predictions. However there is one more factor to the new potential, which also makes chemical sense: the parameters are based on all bonds to hydrogen atoms being shortened by 0.1 Å. Although the nuclear position may be well determined ab initio or by neutron diffraction, the electrons, that is, the point at which the forces act, are located between the hydrogen nucleus and the bonded atom.
To cope with these new types there have been several changes to the now familiar Neighbours and DMAREL codes. Usefully, but not most importantly to the user, NEIGHCRYS will determine the atomic type and use the appropriate atom-atom potential. Also, if required, it will foreshorten the hydrogen bond lengths by 0.1 Å.

Now, you may be thinking that there is at least one problem caused by this, and you would be correct. Firstly, doesn’t the centre of mass change if you adjust the position of the hydrogen atoms? And aren’t the atomic positions relative to the centre of mass? How can this be right?
Well, this is true, and it has an important implication for our work. An ab initio package, like GAUSSIAN [Gaussian 03, Frisch et al, Gaissian Inc.], needs to know the position of the nuclei because it will deal with the electrons itself. But, and very importantly, we always use NEIGHCRYS with our defined axis system to give us the atomic coordinates relative to the centre of mass, so that our input structure from the SHELX res file, our axes and our multipole moments are all consistent. Thankfully, NEIGHCRYS knows what we need, and the coordinates of the atomic nuclei are given correctly relative to the adjusted centre of mass.

4. How to do a proper run using the WILLIAMS potential

This requires some modified versions of the programs you know. All of the programs you need for this Williams run are supplied in the download you will have received. Below are the steps for minimizing a single crystal structure.

Gather your SHELX file, your new bondlengths file, and your Williams potential file (will01.pots)
Run NEIGHCRYS and answer the questions as normal, with the files that you have. NB: There are two new questions. You DO want to use the Williams potential, and you DO want to use foreshortened hydrogen positions.
Create your dmarel.axis file using the atom labels in fort.21
Run NEIGHCRYS with your new axis file
Check that all your atoms are correctly typed in fort.21. Copy the coordinates for molecule 1 in Å only, as you normally would, and put them in your geom file. Do not be phased by the extra section for hydrogens.

Local axis set for molecule: 1 x=> 0.356964 0.654013 0.666966 y=> 0.448954 -0.746252 0.491476 z=> 0.819156 0.123998 -0.560007

Hydrogen positions have been foreshortened. Positions of all atoms (with hydrogens not foreshortened) in the local axis system and centre of mass of the molecule with foreshortened hydrogens for molecule 1
basis No. Species x y z (Angstroms) Mass
1 C_W3_1____ -0.766275 -0.000279 0.000774 12.010700 3 C_W3_2____ 0.766774 -0.000279 0.000774 12.010700 5 O_W1_1____ -1.433231 1.040756 -0.005307 15.999400 7 O_W1_2____ 1.433269 -1.041596 0.000774 15.999400 9 N_W3_1____ -2.619735 -1.511298 0.112844 14.006740 11 N_W3_2____ -1.231180 -1.259466 0.106814 14.006740 13 N_W3_3____ 1.231738 1.259770 -0.095336 14.006740 15 N_W3_4____ 2.619737 1.511681 -0.119590 14.006740 17 H_W4_1____ -2.743490 -2.486029 -0.134621 1.007940 19 H_W4_2____ -2.944737 -1.408808 1.070991 1.007940 21 H_W4_3____ -0.550925 -2.005695 0.052412 1.007940 23 H_W4_4____ 0.548220 2.003881 -0.052643 1.007940 25 H_W4_5____ 2.928963 1.423745 -1.084509 1.007940 27 H_W4_6____ 2.746009 2.484353 0.134988 1.007940
basis No. Species x y z (AU) Mass
1 C_W3_1____ -1.448050 -0.000528 0.001463 12.010700 3 C_W3_2____ 1.448993 -0.000528 0.001463 12.010700 5 O_W1_1____ -2.708415 1.966745 -0.010029 15.999400 7 O_W1_2____ 2.708487 -1.968333 0.001463 15.999400 9 N_W3_1____ -4.950583 -2.855940 0.213244 14.006740 11 N_W3_2____ -2.326594 -2.380046 0.201850 14.006740 13 N_W3_3____ 2.327649 2.380621 -0.180159 14.006740 15 N_W3_4____ 4.950586 2.856664 -0.225992 14.006740 17 H_W4_1____ -5.184447 -4.697915 -0.254397 1.007940 19 H_W4_2____ -5.564748 -2.662263 2.023881 1.007940 21 H_W4_3____ -1.041098 -3.790215 0.099044 1.007940 23 H_W4_4____ 1.035986 3.786787 -0.099481 1.007940 25 H_W4_5____ 5.534941 2.690489 -2.049426 1.007940 27 H_W4_6____ 5.189208 4.694749 0.255091 1.007940
foreshortened hydrogen atom positions in the same local axis system basis No. Species x y z (Angstroms) Mass
17 H_W4_1____ -2.731276 -2.389829 -0.110198 1.007940 19 H_W4_2____ -2.912778 -1.418886 0.976773 1.007940 21 H_W4_3____ -0.618196 -1.931900 0.057792 1.007940 23 H_W4_4____ 0.615808 1.930301 -0.056865 1.007940 25 H_W4_5____ 2.898560 1.432391 -0.989636 1.007940 27 H_W4_6____ 2.733548 2.388366 0.109865 1.007940
basis No. Species x y z (AU) Mass
17 H_W4_1____ -5.161366 -4.516125 -0.208244 1.007940 19 H_W4_2____ -5.504355 -2.681308 1.845834 1.007940 21 H_W4_3____ -1.168221 -3.650763 0.109211 1.007940 23 H_W4_4____ 1.163709 3.647741 -0.107458 1.007940 25 H_W4_5____ 5.477486 2.706827 -1.870143 1.007940 27 H_W4_6____ 5.165659 4.513359 0.207615 1.007940>

Notice that I have used this time the Williams potential. These example labels differ from those in the previous section, allowing you to see the different types of labels. Obviously, if there were no foreshortening of the bonds, then the last two sections would be omitted. Please read the text in this example and realise what these coordinates are. You may have realised that the top block of atomic positions is what must be copied and placed in your mol.geom file, and used for Gaussian. The lower block shows the hydrogen nuclei at their foreshortened position, as used by DMACRYS with the Williams potential.
Use these same atomic positions to generate your .com file for Gaussian, and run the job
You need a simple GDMA input file as before. You may use the same input file as you have used before. However, you must use not gdma but gdma_for. This will calculate the multipoles at the shortened hydrogen positions, using GDMA2.2 (available from Prof. Anthony Stone's group).
At this point you would usually make a cadpac.charges file using gdmaneigh. Instead, use gdmaneighcrys (supplied with your download). Your answers will be as usual. This code will test the atomic coordinates of all atoms between your geom file and your DMA file. They must be the same (they should be unless you have made an error), but it allows the hydrogen positions to be exactly 0.1 Å shorter.
Now you have everything you need. Run NEIGHCRYS again, remembering to answer that you are using the Williams potential with foreshortening of bonds to hydrogen atoms. This will give you’re your usual dmain file.

5. How to do a ‘wrong’ Williams run

You may have realised that it is possible to use the Williams potential without foreshortening the hydrogen positions. This is an incorrect use of the Williams potential. It should be clear how to do this. Please don’t, unless you have a good reason.

Enter secure pages (For project members only - password required)

Warning: This page uses Cascading Style Sheets. We recommend you upgrade to a current version of your favourite browser.