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.mol2This 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):
Wit!P> mpeoe
MPEOE> edit
MPEOE resetThe 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
MPEOE atom C {10*,1,1}
MPEOE atom O {10*,1,1}
MPEOE atom H {10,1,1}
MPEOE bond RR 0.5*
MPEOE> tolerance 0.0000001This 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> 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*
MPEOE> maxF 0.2The show errors command may now be used to find out which net atomic charges are not reproduced well by the current parameter set:
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*
MPEOE> show -print error 0.06 AM1The 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):
/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
MPEOE atom CCO {10.335*,1.000,1.000}and do a few more cycles of refinement:
MPEOE> refine AM1This 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.
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