MPEOE (parameter fitting, exercise 2)


In this exercise, we will develop a simple set of MPEOE from scratch. The file am1.mol2 contains a set of 61 molecules with conformations and net atomic charges obtained with MOPAC/AM1 (full geometry optimization). The molecules contain only C, H, and O atoms, and all bonds are single bonds.

Start by reading the molecules with the target charges, and replace the current MPEOE parameters by a very simple parameter set:

Wit!P> read mol2 am1.mol2
Wit!P> mpeoe
MPEOE> edit
This will call the vi-editor (or the editor ot your choice, as defined by the WNP_EDIT environment variable). Replace the current parameters by the following lines (it is assumed that you are familiar with the swiles notation used in MPEOE parameter definitions):
MPEOE reset
MPEOE atom C {10*,1,1}
MPEOE atom O {10*,1,1}
MPEOE atom H {10,1,1}
MPEOE bond RR 0.5*
The default flavor of swiles in this context is hard. Write the file back to disk, and exit from the editor. The parameters will automatically be read back into Wit!P. The current model is extremely simple, containing just three variable parameters (marked *). The  parameter of H has been kept fixed, since the  parameters are determined only up to an additive constant. All net atomic charges computed with this parameter set are zero. If you don't believe this, you may use the show qcalc or show -print qcalc command to view the charges calculated with the current parameter set. Let's do a few cycles of parameter refinement:
MPEOE>  tolerance 0.0000001
MPEOE>  report 100
MPEOE>  refine AM1
Number of atoms: 1173
MPEOE refinement using 3 parameters:
iter     0  0.12740  0.18521 0.0578038
iter    10  0.03642  0.05279 0.0163716
iter    20  0.03282  0.03306 0.0002432
iter    30  0.03268  0.03269 0.0000107
iter    40  0.03267  0.03268 0.0000074
iter    50  0.03262  0.03264 0.0000127
iter    60  0.03261  0.03261 0.0000006
rms=0.032613
MPEOE> list
MPEOE atom C {10.224*,1.000,1.000}
MPEOE atom O {10.461*,1.000,1.000}
MPEOE bond RR 0.587*
This is a reasonable start, given that we have only two atom and one bond parameter in the refinement, but it is not necessarily the best three parameter model. Indeed the following trick will yield a lower rms after an additional few refinement cycles:
MPEOE>  maxF 0.2
MPEOE>  refine AM1
Number of atoms: 1173
MPEOE refinement using 3 parameters:
iter     0  0.05695  0.13729 0.0803399
iter    10  0.01982  0.02723 0.0074130
iter    20  0.01535  0.01554 0.0001888
iter    30  0.01526  0.01527 0.0000060
rms=0.015259
MPEOE> maxF 0.85
MPEOE> refine AM1
Number of atoms: 1173
MPEOE refinement using 3 parameters:
iter     0  0.01526  0.16037 0.1451105
iter    10  0.01526  0.01608 0.0008231
iter    20  0.01520  0.01523 0.0000260
iter    30  0.01520  0.01520 0.0000007
rms=0.015199
MPEOE> list
MPEOE atom C {10.335*,1.000,1.000}
MPEOE atom O {10.850*,1.000,1.000}
MPEOE bond RR 0.218*
The show errors command may now be used to find out which net atomic charges are not reproduced well by the current parameter set:
MPEOE>  show -print error 0.06 AM1
/AM1/037/C6       -0.147 = -0.242 +0.094
/AM1/040/C5       -0.143 = -0.234 +0.091
/AM1/040/C13      -0.143 = -0.231 +0.088
/AM1/044/C2        0.025 = -0.046 +0.070
/AM1/044/C4        0.025 = -0.048 +0.073
/AM1/056/C3        0.023 = -0.057 +0.080
/AM1/057/C4       -0.143 = -0.216 +0.073
/AM1/058/C4       -0.143 = -0.215 +0.072
8 atoms have abs(dq)>0.0600
The 8 atoms reported by the show errors command share a common feature: they are all described by the swiles CCO. It is therefore reasonable to attempt to improve the model by the introduction of extra parameters for this fragment. Use the edit command to insert the following line into the current parameter set (before the  MPEOE atom C  definition):
MPEOE atom CCO {10.335*,1.000,1.000}
and do a few more cycles of refinement:
MPEOE> refine AM1
Number of atoms: 1173
MPEOE refinement using 4 parameters:
iter     0  0.01520  0.18028 0.1650800
iter    10  0.01508  0.01976 0.0046843
iter    20  0.01088  0.01189 0.0010051
iter    30  0.01076  0.01081 0.0000501
iter    40  0.01073  0.01074 0.0000038
iter    50  0.01073  0.01073 0.0000003
iter    60  0.01073  0.01073 0.0000003
rms=0.010731
MPEOE> list
MPEOE atom CCO {10.357*,1.000,1.000}
MPEOE atom C {10.303*,1.000,1.000}
MPEOE atom O {10.788*,1.000,1.000}
MPEOE bond RR 0.237*
MPEOE> lab rem *
MPEOE> show -print err 0.06 *
0 atoms have abs(dq)>0.0600
This has resulted in a four parameter model, with a root mean square difference between MOPAC/AM1 and MPEOE charges of 0.011 and no errors larger than 0.06. Of course, this model can easily be improved by the introduction of additional parameters. Once a satisfactory model is obtained, it may be written to a disk file with the MPEOE write command.



A.Widmer, NIBR/CPC/CSG-SB