Re: AMBER: An issue with the "multimolecule and multiconformational" resp fit

From: FyD <fyd.u-picardie.fr>
Date: Tue, 11 Apr 2006 21:39:40 +0200

Quoting Jason K <jmk9jmk.gmail.com>:

> I am in the process of building amino-acid variants for simulations. Our
> interest lies in protonated/deprotonated versions of existing amino acids
> for "gas-phase" calculations. After some testing on an arginine residue, we
> decided that calculations will be run at a B3LYP/6-31G* level. In any other
> respect I am trying to follow the J. Comp. Chem. 1995 paper as well as
> possible.

B3LYP/6-31G* for geometry optimization or MEP computation ?

> The first problem that we have encountered was building the starting
> geometries stated in table II, p. 1362 in the aforementioned paper. It would
> really help if the coordinated for these geometries were available.
> Furthermore, how significant is it to start the calculations from such a
> geometry? Especially for the 'idealised' alpha-helical torsions (-60, -40),
> the conformation of backbone atoms is not retained after optimisation. The
> C5 geometry seems to 'survive' more frequently.

'Often' two conformations (one close to alpha helix and the second to extended
forms) are used in a RESP fit. To get a good starting structure for the target
conformation a good way is to first use (phi, psi) dihedral constaints in the
geom. opt. step. Then, reoptimize the structures removing these constraints and
see how different is the structure... It depends on the aminoacids studied.
Another problem is with long side chains with many different chi dihedrals and
potential hydrogen bonds. Sometimes you cannot avoid the use of dihedral
constraints during the geom. opt. step (because the conformation without these
constraints is not 'alpha-helix like' or 'extended like' or because you get
H-bond without them).

Good examples were just submitted by J. Pecher in R.E.DD.B. about Central
fragments of Cysteine derivatives

1) Cys-SO3-: Geom. opt. with (phi, psi) dihedral constraints (to avoid HBond)
 used during the geom. opt. step and using a _single_ conformation RESP fit (the
other conf. had one HBond)
=> See the F-65 projects for Cornell et al. FF, and the F-66 project for the
Duan et al. force field

2) Cys-SSO3-: Geom. opt. with (phi, psi) dihedral constraints (to avoid HBond)
used during the geom. opt. step and using _two_ conformation RESP fit
=> See the F-67 projects for Cornell et al. FF and the F-68 project for the Duan
et al. force field

3) Cys-SMe: Geom. opt. without (phi, psi) dihedral constraints used during the
geom. opt. step and using _two_ conformation RESP fit.
=> See the F-69 projects for Cornell et al. FF and the F-70 project for the Duan
et al. force field

See the abstract in each project for more information about the
orientation/conformation/constraints selected in the charge derivation procedure
http://www.u-picardie.fr/labo/lbpd/REDDB/up/F-65/ up to
http://www.u-picardie.fr/labo/lbpd/REDDB/up/F-70/

Please see also the new "List project" tool of R.E.DD.B. just released:
Project-code Project-name Author Submission date
F-70
Central fragment of S-methyl-L-cysteine Julien Pecher 11-04-2006 13:17
F-69
Central fragment of S-methyl-L-cysteine Julien Pecher 11-04-2006 13:13
F-68
Central fragment of S-sulfonate-L-cysteine Julien Pecher 11-04-2006 12:51
F-67
Central fragment of S-sulfonate-L-cysteine Julien Pecher 11-04-2006 12:39
F-66
Central fragment of sulfonate-L-cysteine Julien Pecher 11-04-2006 10:47
F-65
Central fragment of sulfonate-L-cysteine Julien Pecher 11-04-2006 10:37

> Now, my main question is about the multi-molecular RESP fitting method. I
> have followed all instructions available to extract the ESP data for two
> conformations of the residue (it is a C-terminal amidated Alanine) from the
> Gaussian output files, catenated the two, then generated resp input files
> for the fit. The input for the first fitting step looks like this:

Please look for C-terminal residues in the "download section" of R.E.DD.B.
you will find two examples : alpha-aminoisobutyric acid and O-Methyl-tyrosine

Using "Result(s) for search by Molecule keyword C-terminal"
You will get many examples of RESP inputs for C-terminal AA

Using "Result(s) for search by Molecule keyword N-terminal"
You will get many examples of RESP inputs for N-terminal AA

You can now upload a single .tar.bz2 file by project instead of all the files
one after the others...

> Resp charges for organic molecule
>
> &cntrl
>
> nmol = 2,
> ihfree = 1,
> ioutopt = 1,
>
> &end
???????????????
> 130.0
??????????????????
See
http://www.u-picardie.fr/labo/lbpd/REDDB/uploadfile/F-11/input1.in for instance
> Resp charges for organic molecule
> 0 19
> 1 0
> 6 0
> 1 0
> 1 0
> 6 0
> 8 0
> 7 0
> 1 0
> 6 0
> 1 0
> 6 0
> 8 0
> 6 0
> 1 0
> 1 0
> 1 0
> 7 0
> 1 0
> 1 18
>
> 130.0
See
http://www.u-picardie.fr/labo/lbpd/REDDB/uploadfile/F-11/input1.in for instance
> Resp charges for organic molecule
> 0 19
> 1 0
> 6 0
> 1 0
> 1 0
> 6 0
> 8 0
> 7 0
> 1 0
> 6 0
> 1 0
> 6 0
> 8 0
> 6 0
> 1 0
> 1 0
> 1 0
> 7 0
> 1 0
> 1 18

Why two lines?

>
> 6 0
> 1 1 1 2 1 3 1 4 1 5 1 6

Jump a line here See
http://www.u-picardie.fr/labo/lbpd/REDDB/uploadfile/F-11/input1.in for instance

> 2
> 1 1 2 1
> 2
> 1 2 2 2
> 2
> 1 3 2 3
> 2
> 1 4 2 4
> 2
> 1 5 2 5
> 2
> 1 6 2 6
> 2
> 1 7 2 7
> 2
> 1 8 2 8
> 2
> 1 9 2 9
> 2
> 1 10 2 10
> 2
> 1 11 2 11
> 2
> 1 12 2 12
> 2
> 1 13 2 13
> 2
> 1 14 2 14
> 2
> 1 15 2 15
> 2
> 1 16 2 16
> 2
> 1 17 2 17
> 2
> 1 18 2 18
> 2
> 1 19 2 19
>
>
>
>
>
>
> With the constraints I'm trying to achieve the following:
> - Make all 'identical' atoms from the two conformations equivalent
> - Impose a net charge of zero to the capping group (an acetyl moiety - first
> six atoms)
>
> Unfortunately though, with the current setup I get the following output
> file:
>
> -----------------------------------------------
> Restrained ESP Fit 2.3 Amber 4.1
> -----------------------------------------------
> Resp charges for organic molecule
> -----------------------------------------------
>
>
> inopt = 0 ioutopt = 1
> nmol = 2 iqopt = 0
> ihfree = 1 irstrnt = 1
> iunits = 0 qwt = 0.00050000
>
> %RESP-I-MULT_MOL, multiple-molecule run of 2 molecules
>
> Reading input for molecule 1 weight: 130.000
> Resp charges for organic molecule
>
> Total charge (ich): 0
> Number of centers: 19
> 1 1 0
> 2 6 0
> 3 1 0
> 4 1 0
> 5 6 0
> 6 8 0
> 7 7 0
> 8 1 0
> 9 6 0
> 10 1 0
> 11 6 0
> 12 8 0
> 13 6 0
> 14 1 0
> 15 1 0
> 16 1 0
> 17 7 0
> 18 1 0
> 19 1 18
>
>
> Reading input for molecule 2 weight: 130.000
> Resp charges for organic molecule
>
> Total charge (ich): 0
> Number of centers: 19
> 20 1 0
> 21 6 0
> 22 1 0
> 23 1 0
> 24 6 0
> 25 8 0
> 26 7 0
> 27 1 0
> 28 6 0
> 29 1 0
> 30 6 0
> 31 8 0
> 32 6 0
> 33 1 0
> 34 1 0
> 35 1 0
> 36 7 0
> 37 1 0
> 38 1 18
> --------------------------------
> reading mult_mol constraint info
> --------------------------------
> 1 1 1 2 1 3 1 4 1 5 1 6
> 1 1 2 1
> 1 2 2 2
> 1 3 2 3
> 1 4 2 4
> 1 5 2 5
> 1 6 2 6
> 1 7 2 7
> 1 8 2 8
> 1 9 2 9
> 1 10 2 10
> 1 11 2 11
> 1 12 2 12
> 1 13 2 13
> 1 14 2 14
> 1 15 2 15
> 1 16 2 16
> 1 17 2 17
> 1 18 2 18
> 1 19 2 19
>
> --------------------
> Atom Ivary
> --------------------
> 1 0
> 6 1
> 1 1
> 1 1
> 6 1
> 8 1
> 7 0
> 1 0
> 6 0
> 1 0
> 6 0
> 8 0
> 6 0
> 1 0
> 1 0
> 1 0
> 7 0
> 1 0
> 1 18
>
> 1 1
> 6 1
> 1 1
> 1 1
> 6 1
> 8 1
> 7 7
> 1 8
> 6 9
> 1 10
> 6 11
> 8 12
> 6 13
> 1 14
> 1 15
> 1 16
> 7 17
> 1 18
> 1 18
>
>
> ----------------------------------------------------------------------------
>
>
> Total number of atoms = 38
> Weight factor on initial charge restraints= 0.000500
>
>
> There are 2 charge constraints
>
> Reading esp"s for molecule 1
> total number of atoms = 19
> total number of esp points = 5800
>
>
> center X Y Z
> 1 0.5244048E+01 -0.1457098E+01 -0.3127389E+01
> 2 0.5739154E+01 -0.7333327E+00 -0.1252781E+01
> 3 0.6891891E+01 -0.2131726E+01 -0.2625642E+00
> 4 0.6899441E+01 0.9674244E+00 -0.1487107E+01
> 5 0.3460143E+01 -0.8705212E-01 0.3761632E+00
> 6 0.3643445E+01 0.2927828E+00 0.2647614E+01
> 7 0.1184913E+01 0.1687525E-01 -0.8880635E+00
> 8 0.1222708E+01 -0.4359598E-01 -0.2798576E+01
> 9 -0.1099770E+01 0.1024093E+01 0.2911255E+00
> 10 -0.6254483E+00 0.1184722E+01 0.2301794E+01
> 11 -0.3359876E+01 -0.8108380E+00 0.1229399E+00
> 12 -0.5531173E+01 -0.7007104E-01 0.3780529E+00
> 13 -0.1812204E+01 0.3633803E+01 -0.7387752E+00
> 14 -0.3565871E+01 0.4272525E+01 0.1399474E+00
> 15 -0.3060959E+00 0.4994410E+01 -0.3381532E+00
> 16 -0.2120229E+01 0.3569551E+01 -0.2787238E+01
> 17 -0.2757046E+01 -0.3282598E+01 -0.3022484E+00
> 18 -0.9447968E+00 -0.3866517E+01 -0.2002032E+00
> 19 -0.4153550E+01 -0.4571395E+01 -0.1057169E+00
>
> Reading esp"s for molecule 2
> total number of atoms = 19
> total number of esp points = 5736
>
>
> center X Y Z
> 1 -0.4955020E+01 -0.3620564E+01 0.1067614E+01
> 2 -0.5644746E+01 -0.1866889E+01 0.2134578E+00
> 3 -0.6614181E+01 -0.2314741E+01 -0.1557216E+01
> 4 -0.7037462E+01 -0.1005155E+01 0.1477685E+01
> 5 -0.3567910E+01 0.4737354E-01 -0.3232244E+00
> 6 -0.4013856E+01 0.2116631E+01 -0.1277536E+01
> 7 -0.1188755E+01 -0.6802031E+00 0.3362900E+00
> 8 -0.7938266E+00 -0.2445213E+01 0.9674584E+00
> 9 0.9995687E+00 0.8787887E+00 -0.1909436E+00
> 10 0.7558073E+00 0.1729169E+01 -0.2067442E+01
> 11 0.3308789E+01 -0.8805783E+00 -0.2041717E+00
> 12 0.3218052E+01 -0.3031085E+01 0.6537639E+00
> 13 0.1311404E+01 0.3025514E+01 0.1755474E+01
> 14 0.2929026E+01 0.4233026E+01 0.1300050E+01
> 15 -0.3987813E+00 0.4182050E+01 0.1740356E+01
> 16 0.1587294E+01 0.2246943E+01 0.3652759E+01
> 17 0.5463092E+01 0.1360640E+00 -0.1171711E+01
> 18 0.5499021E+01 0.1844376E+01 -0.2014529E+01
> 19 0.7031549E+01 -0.9505511E+00 -0.1232183E+01
> Initial ssvpot = 74200.412
>
>
> Number of unique UNfrozen centers= 13
>
>
> Non-linear optimization requested.
> qchnge = 0.8392109792E-01
> qchnge = 0.5704251783E-07
> qchnge = 0.3061276458E-12
>
> Convergence in 2 iterations
>
>
>
> Resp charges for organic molecule
>
> Point Charges Before & After Optimization
>
> no. At.no. q(init) q(opt) ivary d(rstr)/dq
> 1 1 0.000000 0.001696 0 0.000000
> 2 6 0.000000 0.001696 1 0.004999
> 3 1 0.000000 0.001696 1 0.000000
> 4 1 0.000000 0.001696 1 0.000000
> 5 6 0.000000 0.001696 1 0.004999
> 6 8 0.000000 0.001696 1 0.004999
> 7 7 0.000000 -0.536747 0 0.000916
> 8 1 0.000000 0.507795 0 0.000000
> 9 6 0.000000 -0.297137 0 0.001595
> 10 1 0.000000 -0.036016 0 0.000000
> 11 6 0.000000 1.199541 0 0.000415
> 12 8 0.000000 -0.669359 0 0.000739
> 13 6 0.000000 0.245507 0 0.001886
> 14 1 0.000000 0.000739 0 0.000000
> 15 1 0.000000 -0.049329 0 0.000000
> 16 1 0.000000 -0.059168 0 0.000000
> 17 7 0.000000 -1.386707 0 0.000360
> 18 1 0.000000 0.535351 0 0.000000
> 19 1 0.000000 0.535351 18 0.000000
>
> 20 1 0.000000 0.001696 1 0.000000
> 21 6 0.000000 0.001696 1 0.004999
> 22 1 0.000000 0.001696 1 0.000000
> 23 1 0.000000 0.001696 1 0.000000
> 24 6 0.000000 0.001696 1 0.004999
> 25 8 0.000000 0.001696 1 0.004999
> 26 7 0.000000 -0.536747 7 0.000916
> 27 1 0.000000 0.507795 8 0.000000
> 28 6 0.000000 -0.297137 9 0.001595
> 29 1 0.000000 -0.036016 10 0.000000
> 30 6 0.000000 1.199541 11 0.000415
> 31 8 0.000000 -0.669359 12 0.000739
> 32 6 0.000000 0.245507 13 0.001886
> 33 1 0.000000 0.000739 14 0.000000
> 34 1 0.000000 -0.049329 15 0.000000
> 35 1 0.000000 -0.059168 16 0.000000
> 36 7 0.000000 -1.386707 17 0.000360
> 37 1 0.000000 0.535351 18 0.000000
> 38 1 0.000000 0.535351 18 0.000000
>
> Sum over the calculated charges: 0.000
>
> Statistics of the fitting:
> The initial sum of squares (ssvpot) 74200.412
> The residual sum of squares (chipot) 74197.016
> The std err of estimate (sqrt(chipot/N)) 3.59657
> ESP relative RMS (SQRT(chipot/ssvpot)) 0.99998
>
>
> Obviously, the charge constraint for the first 6 atoms does not set the net
> charge to zero, but 'forces' every atom to have a charge of (almost) zero.

Yes because you did not jump a line at the right position...
See
http://www.u-picardie.fr/labo/lbpd/REDDB/uploadfile/F-11/input1.in for instance

> This makes the fit a bit "irregular" (attributing e.g. unusually high
> charges to particular atoms). Is there a way to correct this behaviour? I
> have tried different 'spacing' versions of my input file but this is the
> best result I've got so far, which is obviously not good enough. I will be
> very grateful if anyone could help me on this matter!
> I will attached the esp file for your reference. Any advice will be greatly
> appreciated.

You have all the information you need in R.E.DD.B.

regards, Francois

-----------------------------------------------------------------------
http://www.rosswalker.co.uk/movie_sounds/
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Apr 12 2006 - 06:07:11 PDT
Custom Search