From rabe@fu-berlin.de Fri Nov 1 06:18:00 1996 Received: from ki1.chemie.fu-berlin.de for rabe@fu-berlin.de by www.ccl.net (8.8.2/950822.1) id FAA20987; Fri, 1 Nov 1996 05:46:23 -0500 (EST) Received: by ki1.chemie.fu-berlin.de (Smail3.1.28.1) from Bose.Chemie.FU-Berlin.DE (160.45.26.10) with smtp id ; Fri, 1 Nov 96 11:46 MET Received: by Bose.Chemie.FU-Berlin.DE (Smail3.2) from Bose.Chemie.FU-Berlin.DE (127.0.0.1) with smtp id ; Fri, 1 Nov 1996 11:46:17 +0100 (MET) Sender: rabe@www.ccl.net Message-ID: <3279D4F9.15FB@fu-berlin.de> Date: Fri, 01 Nov 1996 11:46:17 +0100 From: Bjoern Rabenstein Organization: Freie Universitaet Berlin, Germany X-Mailer: Mozilla 2.01S (X11; I; IRIX64 6.2 IP26) MIME-Version: 1.0 To: Computational Chemistry Mailing List Subject: Fitting point charges to ESP with Spartan Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Hi! Does anybody know the exact method used by Spartan to calculate point charges by fitting to the electrostatic potential? Or at least a reference to relevant papers? Nothing is written about that in the manual (except some very general information about the fitting procedure). The program tells only the number of points used. I did a calculation with Gaussian using CHELPG and got similar, but not identical results as with Spartan. In a book about using Spartan I found a reference to the Breneman & Wiberg-paper (Journal of Comp. Chem, 11, pp. 361-373 (1990), but I'm not very sure that Spartan is using exact the method described in the paper. Thanks for help! -- Bjoern Rabenstein * Freie Universitaet Berlin * Inst. f. Biochemie Inst. f. Kristallographie * AG Knapp * Takustrasse 6 * D-14195 Berlin [Telephon] +49-30-838-3484 [Telefax] +49-30-838-3464 [email] rabe@fu-berlin.de [WWW] http://www.chemie.fu-berlin.de/~rabe/ From smb@smb.chem.niu.edu Fri Nov 1 13:18:09 1996 Received: from mp.cs.niu.edu for smb@smb.chem.niu.edu by www.ccl.net (8.8.2/950822.1) id MAA22848; Fri, 1 Nov 1996 12:21:30 -0500 (EST) Received: from cz2.chem.niu.edu by mp.cs.niu.edu with SMTP id AA28192 (5.67b/IDA-1.5 for <@mp.cs.niu.edu.chem.niu.edu:CHEMISTRY@www.ccl.net>); Fri, 1 Nov 1996 11:21:31 -0600 Received: from cz.chem.niu.edu by cz2.chem.niu.edu via SMTP (920330.SGI/890607.SGI) (for @mp.cs.niu.edu.chem.niu.edu:CHEMISTRY@www.ccl.net) id AA05625; Fri, 1 Nov 96 08:14:21 -0600 Received: from smb.chem.niu.edu by cz.chem.niu.edu via SMTP (920330.SGI/890607.SGI) (for @cz2.chem.niu.edu:CHEMISTRY@www.ccl.net) id AA25887; Fri, 1 Nov 96 08:14:19 -0600 Received: by smb.chem.niu.edu (940816.SGI.8.6.9/940406.SGI) for CHEMISTRY@www.ccl.net id IAA02588; Fri, 1 Nov 1996 08:14:17 -0600 Date: Fri, 1 Nov 1996 08:14:17 -0600 From: smb@smb.chem.niu.edu (Steven Bachrach) Message-Id: <199611011414.IAA02588@smb.chem.niu.edu> To: CHEMISTRY@www.ccl.net Subject: ECCC3 is now open One last announcement concerning the 3rd Electronic Computational Chemistry Conference (ECCC-3) The conference is now open to all (at no cost). Access to the conference is at http://hackberry.chem.niu.edu/ECCC3 From here you can get to the papers and the discussion sections. A European mirror site is also available at http://rchs1.chemie.uni-regensburg.de:8000/ and I would encourage all European participants to use this site. Welcome to all! Steve Steven Bachrach Department of Chemistry Northern Illinois University DeKalb, Il 60115 Phone: (815)753-6863 smb@smb.chem.niu.edu Fax: (815)753-4802 From huang@nissan.wavefun.com Fri Nov 1 18:18:06 1996 Received: from volvo.wavefun.com for huang@nissan.wavefun.com by www.ccl.net (8.8.2/950822.1) id RAA24566; Fri, 1 Nov 1996 17:47:17 -0500 (EST) Received: by volvo.wavefun.com (931110.SGI/930416.SGI) for chemistry@www.ccl.net id AA29937; Fri, 1 Nov 96 14:48:03 -0800 Received: from nissan.wavefun.com(198.147.95.2) by volvo.wavefun.com via smap (V1.3) id sma029935; Fri Nov 1 14:47:59 1996 Received: by nissan.wavefun.com (931110.SGI/940406.SGI.AUTO) for @volvo.wavefun.com:chemistry@www.ccl.net id AA29316; Fri, 1 Nov 96 14:46:11 -0800 From: "Wayne Huang" Message-Id: <9611011446.ZM29314@nissan.wavefun.com> Date: Fri, 1 Nov 1996 14:46:02 -0800 In-Reply-To: Bjoern Rabenstein "CCL:G:Fitting point charges to ESP with Spartan" (Nov 1, 11:46am) References: <3279D4F9.15FB@fu-berlin.de> Reply-To: huang@wavefun.com X-Mailer: Z-Mail (3.2.2 10apr95 MediaMail) To: Bjoern Rabenstein , Computational Chemistry Mailing List Subject: Re: CCL:G:Fitting point charges to ESP with Spartan Cc: sparlist@nissan.wavefun.com Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii On Nov 1, 11:46am, Bjoern Rabenstein wrote: > Subject: CCL:G:Fitting point charges to ESP with Spartan > Hi! > > Does anybody know the exact method used by Spartan to calculate point > charges by fitting to the electrostatic potential? Or at least a > reference to relevant papers? > > Nothing is written about that in the manual (except some very general > information about the fitting procedure). The program tells only the > number of points used. > > I did a calculation with Gaussian using CHELPG and got similar, but not > identical results as with Spartan. > > In a book about using Spartan I found a reference to the Breneman & > Wiberg-paper (Journal of Comp. Chem, 11, pp. 361-373 (1990), but I'm not > very sure that Spartan is using exact the method described in the paper. Spartan does use CHELPG for electrostatic fit charge partition. See: C. Brenneman, K.B. Wiberg, J. Comput. Chem., 11, 361-373 (1990) L.E. Chirlian M.M. Franci, J. Comput. CHem., 8, 894 (1987) As there is discrepancy on resulting charges from G9x and Spartan, that would not be surprising. Although based on the same idea, the implementation could be different. It is somewhat arbitrary in at least two ways, the number of grid points chosen and the how far the region included. In plain English, it basically paints a layer of grease on the vdw surface. Too far beyond the vdw surface would get the charges too small (eventually approach zero). The extent of the exclusion of the inner volume will also affect the outcome (too close to the centers results in charges being too large). I would tend to believe the two implementations should differ by a factor, may or may not be linearly scalable. In the nutshell, regardless which partition scheme you use, stick with it systematically which should give you sensible comparison (hopefully). Have a good one! --Wayne -- +---------------------------------------------------------------+ | Wayne Huang, Ph.D. | 18401 Von Karman, Suite 370 | | Computational Chemist | Irvine, California 92612 | | Wavefunction, Inc. | 714-955-2120 <> 955-2118(fax)| | huang@wavefun.com | Web: http://www.wavefun.com | +---------------------------------------------------------------+ From mn1@helix.nih.gov Fri Nov 1 20:18:08 1996 Received: from helix.nih.gov for mn1@helix.nih.gov by www.ccl.net (8.8.2/950822.1) id TAA24805; Fri, 1 Nov 1996 19:19:44 -0500 (EST) Received: (from mn1@localhost) by helix.nih.gov (8.7.6/8.7.3) id TAA15139 for CHEMISTRY@www.ccl.net; Fri, 1 Nov 1996 19:19:46 -0500 (EST) Date: Fri, 1 Nov 1996 19:19:46 -0500 (EST) From: "M. Nicklaus" Message-Id: <199611020019.TAA15139@helix.nih.gov> To: CHEMISTRY@www.ccl.net Subject: G:SUMMARY: Wavefunction Instability Dear CCL'ers, A few weeks ago I asked the question(s) listed below concerning wavefunction instability encountered in the calculation of a cation complexed with the hydroxyl groups of catechol. I am summarizing the answers I got, slightly edited for the sake of brevity. Thanks to all who responded. The bottom line seems to be: 1. This instability is probably an artifact of the HF method; and 2. since the energy difference is so small anyway, forget about it. I performed the UHF/6-311G** geometry optimization in the meantime, reading in the wavefunction from STABLE=OPT (at the geometry optimized at RHF/6-31G*). The energy change this time, i.e. during geometry optimization, was even smaller (0.25 kcal/mol) than the change when going from SP --> STABLE=OPT (1.45 kcal/mol). The spin contamination increased a bit, from S**2 before annihilation 0.3882, after 0.2298 to S**2 before annihilation 0.4404, after 0.3042 From these results, and the responses, I tend to deduce that this is about as good as it gets (at least at the HF level of theory), and that I can probably reasonably trust the energies I've got so far. If anyone thinks otherwise, please let me know. Marc ------------------------------------------------------------------------ Marc C. Nicklaus Lab. of Medicinal Chemistry e-mail: mn1@helix.nih.gov National Cancer Institute, NIH Phone: (301) 402-3111 Bldg 37, Rm 5B29 Fax: (301) 496-5839 BETHESDA, MD 20892-4255 USA WWW: http://www.nci.nih.gov/intra/lmch/MCNBIO.HTM ------------------------------------------------------------------------ ############################ Original Question ######################### I am calculating, at the HF/6-311G** level in Gaussian 94, the divalent cation Mg++ complexed with the hydroxyl groups of catechol. Both catechol and the ion (should) have a spin multiplicity of 1, so that's what I used for the complex. At the current geometry, which was optimized with a smaller basis set, I performed a wavefunction stability test (keyword "Stable"). I got the message The wavefunction has an RHF -> UHF instability. Running a "Stable=Opt" job produced a stable wavefunction at a slightly lower energy. The differences are not huge: the energy is lowered by 1.45 kcal/mol, and the changes in the charge on the Mg++ ion (both Mulliken and NPA) are less than 10^-3. (Interestingly, the complex of Mg++ with the aromatic ring of the catechol molecule did not show any wavefunction instability.) Using UHF instead of (R)HF, while keeping mult. = 1, produced exactly identical results compared with the HF calculations, both with and without the keyword "Stable=Opt". The spin contamination produced by the wavefunction optimization is Annihilation of the first spin contaminant: S**2 before annihilation 0.3882, after 0.2298 Before I burn weeks' worth of CPU cycles on potentially unnecessary and/or meaningless calculations, I'd like to ask for advice on: (a) Are these results a reason for concern, i.e. does S**2 = 0.3882/0.2298 mean that the optimized wavefunction is useless; or, on the other hand, does the small energy change of only 1.45 kcal/mol mean that I shouldn't even bother optimizing the wavefunction? (b) Since I intend to optimize the geometry of the complex at the HF/6-311G** level, is it possible the the wavefunction instability will go away anyway after that? [I guess I could simply try that, but see (c).] (c) If I (should) use "Stable=Opt", how do I combine that with the geometry optimization? Can you use both in the same job? If so, when does G94 do the wavefunction optimization; i.e., (only) for the initial geometry, and/or for each intermediate geometry, and/or (only) for the final geometry? This would appear to have a major impact on the CPU time used. If you can't combine "Stable=Opt" and geometry optimization in a single job, how should one proceed? (d) Any other comments, ideas, suggestions or warnings concerning this problem? ################################# Responses ############################### From: David Heisterberg Do you mean you did a UHF calculation starting from Gaussian's Huckel guess (or your RHF results)? And did you use GUESS=MIX? Anyway, if you run UHF with the results from the STABLE=OPT run you should recover the lower energy solution. I bring that up because what you might want to do is look at the UHF natural orbitals and I don't know off-hand whether you can get them as a result of an RHF stability calculation. The UHF NOs may tell you more about the nature of the instability. If the interaction between the Mg++ and the OH group is fairly weak at your current geometry the system may not be well described by a single determinant. Look for NOs with populations greater that about 0.02 and less than 0.98. But you may find that the instability goes away as you optimize the geometry with the larger basis. You might kill two birds with one stone by doing a UHF optimization starting from the STABLE output. If the spin contamination gets small at convergence you're probably ok. ------------------------------------------------------------------- From: Frank Jensen Bottom line: don't bother with the UHF wave function. If you want the singlet, run everything RHF. The UHF instability is due to a small amount of electron correlation being included in UHF. If you where to do this at R/U MP2, the RMP2 would be the lowest energy wave function. ------------------------------------------------------------------- From: Ulrike Salzner HF theory has a tendency to produce triplet (and sometimes singlet) instabilities that are artifacts of the method. That can be understood by the neglect of part of the electron correlation in HF theory. Since the Pauli principle is obeyed by using determinatal wave functions, electrons with same spins are correlated. They can not be at the same place. Electrons with different spins, however, are not correlated. That is, the error in the HF approximation is smaller for triplets than for singlets. This can lead to erronous preference for the triplet. Since your system should be a singlet, as you said, and since the energy differences between your singlet and triplet calculations are quite small, I think you can forget about the triplet instability. ------------------------------------------------------------------- From: Doug Fox I suspect that the reason you got back the same energy is because you did not read in the result of the STABILITY calculation as the initial guess. The default initial guess is RHF so just saying UHF is not enough to break the spin symmetry. You might try UHF with GUESS=MIX or read in the solution from STABLE=OPT. Will this instability go away, no option but to check. Normally I would not expect so unless the structure changes dramatically and what you are seeing is due to a bad initial guess at the structure. Since many people do optimizations without ever checking it is not easy to guess how common that is. So I would optimize normally and then run STABLE=OPT again at the end to see if you have stayed on the more stable surface. ------------------------------------------------------------------- From: Michael Braeuer I had to deal with the same problem optimizing a simple problem as acrolein. So far I didn't find any explanation but I'm working on it. ------------------------------------------------------------------- From: Anthony P Scott I would be interested in the general comments on this. I have noticed similar spin contamination values and have not been to sure what to make of them I think that an instability in the wavefunction (given a suitable basis) will not go away with increasing basis set size. I would be happy to be corected on this however., I would proceed (and do) as follows. a. optimize the geometry at rhf or whatever. b. do a stability run with stable=opt. c. with this wavefunction (guess=read) reoptimize the geometry. Generally you will find that the geometry does not change much. ######################## End of Responses #######################