From: <ni...@us...> - 2008-08-27 19:45:48
|
Revision: 9785 http://jmol.svn.sourceforge.net/jmol/?rev=9785&view=rev Author: nicove Date: 2008-08-27 19:45:44 +0000 (Wed, 27 Aug 2008) Log Message: ----------- Code cleanup Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/Minimizer.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java Modified: trunk/Jmol/src/org/jmol/minimize/Minimizer.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2008-08-27 19:43:09 UTC (rev 9784) +++ trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2008-08-27 19:45:44 UTC (rev 9785) @@ -707,7 +707,7 @@ void endMinimization() { updateAtomXYZ(); setMinimizationOn(false); - boolean failed = pFF.DetectExplosion(); + boolean failed = pFF.detectExplosion(); if (failed) restoreCoordinates(); viewer.setIntProperty("_minimizationStep", pFF.getCurrentStep()); Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java 2008-08-27 19:43:09 UTC (rev 9784) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java 2008-08-27 19:45:44 UTC (rev 9785) @@ -410,7 +410,7 @@ atoms[i].coord[j] = coordSaved[i][j]; } - public boolean DetectExplosion() { + public boolean detectExplosion() { for (int i = 0; i < atomCount; i++) { MinAtom atom = atoms[i]; for (int j = 0; j < 3; j++) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2010-04-14 22:51:09
|
Revision: 12861 http://jmol.svn.sourceforge.net/jmol/?rev=12861&view=rev Author: hansonr Date: 2010-04-14 22:51:03 +0000 (Wed, 14 Apr 2010) Log Message: ----------- UFF calc intern Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/Minimizer.java trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java Modified: trunk/Jmol/src/org/jmol/minimize/Minimizer.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2010-04-14 22:34:51 UTC (rev 12860) +++ trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2010-04-14 22:51:03 UTC (rev 12861) @@ -335,7 +335,7 @@ for (int j = 0, pt = 0; j < atomCountFull; j++) if (bsAtoms.get(j)) { if (search.get(j)) { - minAtoms[pt].type = data[1]; + minAtoms[pt].type = data[1].intern(); //System.out.println("pt" +pt + data[1]); } pt++; Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java 2010-04-14 22:34:51 UTC (rev 12860) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java 2010-04-14 22:51:03 UTC (rev 12861) @@ -284,7 +284,7 @@ a = atoms[ia = angle[0]]; b = atoms[ib = angle[1]]; c = atoms[ic = angle[2]]; - boolean isHXH = (a.type.equals("H_") && c.type.equals("H_")); + boolean isHXH = (a.type == "H_" && c.type == "H_"); parA = getParameter(a.type, ffParams); parB = getParameter(b.type, ffParams); @@ -695,7 +695,7 @@ double koop = KCAL6; switch (elemNo) { case 6: // carbon could be a carbonyl, which is considerably stronger - if ((a.type + c.type + d.type).indexOf("O_2") >= 0) { + if (a.type == "O_2" || c.type == "O_2" || d.type == "O_2") { koop += KCAL44; break; }/* else if (b.type.lastIndexOf("R") == 2) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2010-04-27 14:58:09
|
Revision: 12950 http://jmol.svn.sourceforge.net/jmol/?rev=12950&view=rev Author: hansonr Date: 2010-04-27 14:58:00 +0000 (Tue, 27 Apr 2010) Log Message: ----------- version=12.0.RC9_dev # bug fix: UFF calculation fails to recognize C+ and C- Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/Minimizer.java trunk/Jmol/src/org/jmol/minimize/forcefield/UFF.txt Modified: trunk/Jmol/src/org/jmol/minimize/Minimizer.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2010-04-27 13:28:17 UTC (rev 12949) +++ trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2010-04-27 14:58:00 UTC (rev 12950) @@ -332,8 +332,7 @@ else if (search == null) break; else - for (int j = 0, pt = 0; j < atomCountFull; j++) - if (bsAtoms.get(j)) { + for (int j = bsAtoms.nextSetBit(0), pt = 0; j < atomCountFull && j >= 0; j = bsAtoms.nextSetBit(j + 1)) { if (search.get(j)) { minAtoms[pt].type = data[1].intern(); //System.out.println("pt" +pt + data[1]); @@ -532,10 +531,15 @@ search = tokenTypes[TOKEN_ELEMENT_CHARGED]; search[PT_CHARGE].intValue = n; break; + case '-': + search = tokenTypes[TOKEN_ELEMENT_CHARGED]; + search[PT_CHARGE].intValue = -n; + break; } search[PT_ELEMENT].intValue = elemNo; Object v = viewer.evaluateExpression(search); - //System.out.println(smarts + " minimize atoms=" + v.toString()); + if (Logger.debugging) + Logger.debug(smarts + " minimize atoms=" + v.toString()); return (v instanceof BitSet ? (BitSet) v : null); } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/UFF.txt =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/UFF.txt 2010-04-27 13:28:17 UTC (rev 12949) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/UFF.txt 2010-04-27 14:58:00 UTC (rev 12950) @@ -26,7 +26,9 @@ atom [#5] B_2 Trigonal planar boron atom [#5D4] B_3 Tetrahedral boron atom [#6] C_3 Generic sp3 C -atom [C^2] C_2 sp2 non-aromatic C +atom [C^2] C_2 sp2 non-aromatic C= +atom [CD3] C_2 trivalent C (cation) // bh added Jmol 11.12.RC9 +atom [C-1] C_3 trivalent C (anion) // bh added Jmol 11.12.RC9 atom [C^1] C_1 sp hybridized C atom [c] C_R aromatic C atom [#7] N_3 Generic sp3 N @@ -35,7 +37,7 @@ atom [n] N_R aromatic N atom [#8] O_3 generic, sp3 hybridized O #atom [#8] O_3_z sp3 hybridized O for zeolites -atom [O^2] O_2 sp2 hybridized O +atom [O^2] O_2 sp2 hybridized O atom [O^1] O_1 sp hybridized O atom [o] O_R aromatic O atom [#9] F_ generic F This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2012-05-09 18:10:11
|
Revision: 17116 http://jmol.svn.sourceforge.net/jmol/?rev=17116&view=rev Author: hansonr Date: 2012-05-09 18:10:04 +0000 (Wed, 09 May 2012) Log Message: ----------- refactoring Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/Minimizer.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldUFF.java Removed Paths: ------------- trunk/Jmol/src/org/jmol/minimize/forcefield/UFF.txt Modified: trunk/Jmol/src/org/jmol/minimize/Minimizer.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2012-05-09 18:09:32 UTC (rev 17115) +++ trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2012-05-09 18:10:04 UTC (rev 17116) @@ -445,11 +445,6 @@ return pFF; } - public List<String[]> getAtomTypes() { - getForceField(); - return (pFF == null ? null : pFF.getAtomTypes()); - } - /* *************************************************************** * Minimization thead support ****************************************************************/ Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java 2012-05-09 18:09:32 UTC (rev 17115) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java 2012-05-09 18:10:04 UTC (rev 17116) @@ -56,7 +56,7 @@ return calc.getUnit(); } - public abstract List<String[]> getAtomTypes(); + protected abstract List<String[]> getAtomTypes(); protected abstract Map<String, FFParam> getFFParameters(); Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldUFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldUFF.java 2012-05-09 18:09:32 UTC (rev 17115) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldUFF.java 2012-05-09 18:10:04 UTC (rev 17116) @@ -47,6 +47,11 @@ private BitSet bsAromatic; @Override + public void clear() { + bsAromatic = null; + } + + @Override public boolean setModel(Minimizer m, BitSet bsElements, int elemnoMax) { super.setModel(m); Logger.info("minimize: setting atom types..."); @@ -181,7 +186,7 @@ // open UFF.txt URL url = null; - String fileName = "UFF.txt"; + String fileName = "uff/UFF.txt"; BufferedReader br = null; try { if ((url = this.getClass().getResource(fileName)) == null) { @@ -264,7 +269,7 @@ public List<String[]> getAtomTypes() { List<String[]> types = new ArrayList<String[]>(); //!< external atom type rules URL url = null; - String fileName = "UFF.txt"; + String fileName = "uff/UFF.txt"; try { if ((url = this.getClass().getResource(fileName)) == null) { System.err.println("Couldn't find file: " + fileName); @@ -292,7 +297,7 @@ + fileName); } - Logger.info(types.size() + " force field parameters read"); + Logger.info(types.size() + " UFF parameters read"); return (types.size() > 0 ? types : null); } @@ -403,9 +408,4 @@ Token.tokenExpressionEnd}, }; - @Override - public void clear() { - bsAromatic = null; - } - } Deleted: trunk/Jmol/src/org/jmol/minimize/forcefield/UFF.txt =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/UFF.txt 2012-05-09 18:09:32 UTC (rev 17115) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/UFF.txt 2012-05-09 18:10:04 UTC (rev 17116) @@ -1,281 +0,0 @@ -# -# Open Babel file: UFF.prm -# -# Force field parameters for UFF, the Universal Force Field -# Used by OBForceField and OBUFFForceField -# -# J. Am. Chem. Soc. (1992) 114(25) p. 10024-10035. -# The parameters in this file are taken from the UFF implementation in RDKit -# http://rdkit.org/ -# -# Atom typing rules are based on UFF published atom descriptions -# atom [SMARTS] [atomtype] [description] -# These must be sorted according to the element and complexity -# of the type rule (i.e., generic rules first, more specific rules later) -# -# MODIFIED 3/2008 by Bob Hanson for Jmol (bh) -# -# Parameters follow later -# param Atom r1 theta0 x1 D1 zeta Z1 Vi Uj Xi Hard Radius - -atom [#1] H_ Generic hydrogen -atom [#1D2] H_b Bridging hydrogen -atom [#2] He4+4 Helium -atom [#3] Li Lithium -atom [#4] Be3+2 Generic Be -atom [#5] B_2 Trigonal planar boron -atom [#5D4] B_3 Tetrahedral boron -atom [#6] C_3 Generic sp3 C -atom [C^2] C_2 sp2 non-aromatic C= -atom [C+1] C_2+ trivalent C (cation) // bh added Jmol 12.0.RC9 -atom [C-1] C_3 trivalent C (anion) // bh added Jmol 12.0.RC9 -atom [CA1] C_2 allylic C (anion or cation) // bh added Jmol 12.0.RC9 -atom [C^1] C_1 sp hybridized C -atom [c] C_R aromatic C -atom [#7] N_3 Generic sp3 N -atom [NA1] N_2 allylic N or amide // bh added Jmol 12.0.RC9 -atom [N^2] N_2 sp2 non-aromatic N // bh was [ND2], but this improperly treats N-oxides -atom [N^1] N_1 sp hybridized N // bh was [ND1], but this is more specifically sp -atom [n] N_R aromatic N -atom [#8] O_3 generic, sp3 hybridized O -#atom [#8] O_3_z sp3 hybridized O for zeolites -atom [O^2] O_2 sp2 hybridized O -atom [O^1] O_1 sp hybridized O -atom [o] O_R aromatic O -atom [#9] F_ generic F -atom [#10] Ne4+4 -atom [#11] Na -atom [#12] Mg3+2 -atom [#13] Al3 -atom [#14] Si3 -atom [#15+5] P_3+5 formal charge +5 -#atom [#15] P_3+q Organometallic phosphine ligands -atom [#15] P_3+3 generic phosphorus -atom [#16] S_3+2 generic S -atom [#16+4] S_3+4 S+4 -atom [#16+6] S_3+6 S+6 -atom [S^2] S_2 non-aromatic sp2 S -atom [s] S_R aromatic S -atom [#17] Cl -atom [#18] Ar4+4 -atom [#19] K_ -atom [#20] Ca6+2 -atom [#21] Sc3+3 -atom [#22] Ti6+4 generic Ti (6-valent) -atom [#22D3] Ti3+4 -atom [#23] V_3+5 -atom [#24] Cr6+3 -atom [#25] Mn6+2 -atom [#26] Fe6+2 generic Fe (6-valent) -atom [#26D3] Fe3+2 -atom [#27] Co6+3 -atom [#28] Ni4+2 -atom [#29] Cu3+1 -atom [#30] Zn3+2 -atom [#31] Ga3+3 -atom [#32] Ge3 -atom [#33] As3+3 -atom [#34] Se3+2 -atom [#35] Br -atom [#36] Kr4+4 -atom [#37] Rb -atom [#38] Sr6+2 -atom [#39] Y_3+3 -atom [#40] Zr3+4 -atom [#41] Nb3+5 -atom [#42] Mo6+6 generic Mo (6-valent) -atom [#42D3] Mo3+6 trivalent Mo -atom [#43] Tc6+5 -atom [#44] Ru6+2 -atom [#45] Rh6+3 -atom [#46] Pd4+2 -atom [#47] Ag1+1 -atom [#48] Cd3+2 -atom [#49] In3+3 -atom [#50] Sn3 -atom [#51] Sb3+3 -atom [#52] Te3+2 -atom [#53] I_ -atom [#54] Xe4+4 -atom [#55] Cs -atom [#56] Ba6+2 -atom [#57] La3+3 -atom [#58] Ce6+3 -atom [#59] Pr6+3 -atom [#60] Nd6+3 -atom [#61] Pm6+3 -atom [#62] Sm6+3 -atom [#63] Eu6+3 -atom [#64] Gd6+3 -atom [#65] Tb6+3 -atom [#66] Dy6+3 -atom [#67] Ho6+3 -atom [#68] Er6+3 -atom [#69] Tm6+3 -atom [#70] Yb6+3 -atom [#71] Lu6+3 -atom [#72] Hf3+4 -atom [#73] Ta3+5 -atom [#74] W_6+6 generic W (6-valent) -atom [#74D3+4] W_3+4 -atom [#74D3+6] W_3+6 -atom [#75] Re6+5 generic Re (6-valent) -atom [#75D3] Re3+7 trivalent Re -atom [#76] Os6+6 -atom [#77] Ir6+3 -atom [#78] Pt4+2 -atom [#79] Au4+3 -atom [#80] Hg1+2 -atom [#81] Tl3+3 -atom [#82] Pb3 -atom [#83] Bi3+3 -atom [#84] Po3+2 -atom [#85] At -atom [#86] Rn4+4 -atom [#87] Fr -atom [#88] Ra6+2 -atom [#89] Ac6+3 -atom [#90] Th6+4 -atom [#91] Pa6+4 -atom [#92] U_6+4 -atom [#93] Np6+4 -atom [#94] Pu6+4 -atom [#95] Am6+4 -atom [#96] Cm6+3 -atom [#97] Bk6+3 -atom [#98] Cf6+3 -atom [#99] Es6+3 -atom [#100] Fm6+3 -atom [#101] Md6+3 -atom [#102] No6+3 -atom [#103] Lw6+3 - -# Atom r1 theta0 x1 D1 zeta Z1 Vi Uj Xi Hard Radius -param H_ 0.354 180 2.886 0.044 12 0.712 0 0 4.528 6.9452 0.371 -param H_b 0.46 83.5 2.886 0.044 12 0.712 0 0 4.528 6.9452 0.371 -param He4+4 0.849 90 2.362 0.056 15.24 0.098 0 0 9.66 14.92 1.3 -param Li 1.336 180 2.451 0.025 12 1.026 0 2 3.006 2.386 1.557 -param Be3+2 1.074 109.47 2.745 0.085 12 1.565 0 2 4.877 4.443 1.24 -param B_3 0.838 109.47 4.083 0.18 12.052 1.755 0 2 5.11 4.75 0.822 -param B_2 0.828 120 4.083 0.18 12.052 1.755 0 2 5.11 4.75 0.822 -param C_3 0.757 109.47 3.851 0.105 12.73 1.912 2.119 2 5.343 5.063 0.759 -param C_R 0.729 120 3.851 0.105 12.73 1.912 0 2 5.343 5.063 0.759 -param C_2 0.732 120 3.851 0.105 12.73 1.912 0 2 5.343 5.063 0.759 -param C_2+ 0.732 120 6.851 0.105 12.73 1.912 0 2 5.343 5.063 0.759 -param C_1 0.706 180 3.851 0.105 12.73 1.912 0 2 5.343 5.063 0.759 -param N_3 0.7 106.7 3.66 0.069 13.407 2.544 0.45 2 6.899 5.88 0.715 -param N_R 0.699 120 3.66 0.069 13.407 2.544 0 2 6.899 5.88 0.715 -param N_2 0.685 111.2 3.66 0.069 13.407 2.544 0 2 6.899 5.88 0.715 -param N_1 0.656 180 3.66 0.069 13.407 2.544 0 2 6.899 5.88 0.715 -param O_3 0.658 104.51 3.5 0.06 14.085 2.3 0.018 2 8.741 6.682 0.669 -param O_3_z 0.528 146 3.5 0.06 14.085 2.3 0.018 2 8.741 6.682 0.669 -param O_R 0.68 110 3.5 0.06 14.085 2.3 0 2 8.741 6.682 0.669 -param O_2 0.634 120 3.5 0.06 14.085 2.3 0 2 8.741 6.682 0.669 -param O_1 0.639 180 3.5 0.06 14.085 2.3 0 2 8.741 6.682 0.669 -param F_ 0.668 180 3.364 0.05 14.762 1.735 0 2 10.874 7.474 0.706 -param Ne4+4 0.92 90 3.243 0.042 15.44 0.194 0 2 11.04 10.55 1.768 -param Na 1.539 180 2.983 0.03 12 1.081 0 1.25 2.843 2.296 2.085 -param Mg3+2 1.421 109.47 3.021 0.111 12 1.787 0 1.25 3.951 3.693 1.5 -param Al3 1.244 109.47 4.499 0.505 11.278 1.792 0 1.25 4.06 3.59 1.201 -param Si3 1.117 109.47 4.295 0.402 12.175 2.323 1.225 1.25 4.168 3.487 1.176 -param P_3+3 1.101 93.8 4.147 0.305 13.072 2.863 2.4 1.25 5.463 4 1.102 -param P_3+5 1.056 109.47 4.147 0.305 13.072 2.863 2.4 1.25 5.463 4 1.102 -param P_3+q 1.056 109.47 4.147 0.305 13.072 2.863 2.4 1.25 5.463 4 1.102 -param S_3+2 1.064 92.1 4.035 0.274 13.969 2.703 0.484 1.25 6.928 4.486 1.047 -param S_3+4 1.049 103.2 4.035 0.274 13.969 2.703 0.484 1.25 6.928 4.486 1.047 -param S_3+6 1.027 109.47 4.035 0.274 13.969 2.703 0.484 1.25 6.928 4.486 1.047 -param S_R 1.077 92.2 4.035 0.274 13.969 2.703 0 1.25 6.928 4.486 1.047 -param S_2 0.854 120 4.035 0.274 13.969 2.703 0 1.25 6.928 4.486 1.047 -param Cl 1.044 180 3.947 0.227 14.866 2.348 0 1.25 8.564 4.946 0.994 -param Ar4+4 1.032 90 3.868 0.185 15.763 0.3 0 1.25 9.465 6.355 2.108 -param K_ 1.953 180 3.812 0.035 12 1.165 0 0.7 2.421 1.92 2.586 -param Ca6+2 1.761 90 3.399 0.238 12 2.141 0 0.7 3.231 2.88 2 -param Sc3+3 1.513 109.47 3.295 0.019 12 2.592 0 0.7 3.395 3.08 1.75 -param Ti3+4 1.412 109.47 3.175 0.017 12 2.659 0 0.7 3.47 3.38 1.607 -param Ti6+4 1.412 90 3.175 0.017 12 2.659 0 0.7 3.47 3.38 1.607 -param V_3+5 1.402 109.47 3.144 0.016 12 2.679 0 0.7 3.65 3.41 1.47 -param Cr6+3 1.345 90 3.023 0.015 12 2.463 0 0.7 3.415 3.865 1.402 -param Mn6+2 1.382 90 2.961 0.013 12 2.43 0 0.7 3.325 4.105 1.533 -param Fe3+2 1.27 109.47 2.912 0.013 12 2.43 0 0.7 3.76 4.14 1.393 -param Fe6+2 1.335 90 2.912 0.013 12 2.43 0 0.7 3.76 4.14 1.393 -param Co6+3 1.241 90 2.872 0.014 12 2.43 0 0.7 4.105 4.175 1.406 -param Ni4+2 1.164 90 2.834 0.015 12 2.43 0 0.7 4.465 4.205 1.398 -param Cu3+1 1.302 109.47 3.495 0.005 12 1.756 0 0.7 4.2 4.22 1.434 -param Zn3+2 1.193 109.47 2.763 0.124 12 1.308 0 0.7 5.106 4.285 1.4 -param Ga3+3 1.26 109.47 4.383 0.415 11 1.821 0 0.7 3.641 3.16 1.211 -param Ge3 1.197 109.47 4.28 0.379 12 2.789 0.701 0.7 4.051 3.438 1.189 -param As3+3 1.211 92.1 4.23 0.309 13 2.864 1.5 0.7 5.188 3.809 1.204 -param Se3+2 1.19 90.6 4.205 0.291 14 2.764 0.335 0.7 6.428 4.131 1.224 -param Br 1.192 180 4.189 0.251 15 2.519 0 0.7 7.79 4.425 1.141 -param Kr4+4 1.147 90 4.141 0.22 16 0.452 0 0.7 8.505 5.715 2.27 -param Rb 2.26 180 4.114 0.04 12 1.592 0 0.2 2.331 1.846 2.77 -param Sr6+2 2.052 90 3.641 0.235 12 2.449 0 0.2 3.024 2.44 2.415 -param Y_3+3 1.698 109.47 3.345 0.072 12 3.257 0 0.2 3.83 2.81 1.998 -param Zr3+4 1.564 109.47 3.124 0.069 12 3.667 0 0.2 3.4 3.55 1.758 -param Nb3+5 1.473 109.47 3.165 0.059 12 3.618 0 0.2 3.55 3.38 1.603 -param Mo6+6 1.467 90 3.052 0.056 12 3.4 0 0.2 3.465 3.755 1.53 -param Mo3+6 1.484 109.47 3.052 0.056 12 3.4 0 0.2 3.465 3.755 1.53 -param Tc6+5 1.322 90 2.998 0.048 12 3.4 0 0.2 3.29 3.99 1.5 -param Ru6+2 1.478 90 2.963 0.056 12 3.4 0 0.2 3.575 4.015 1.5 -param Rh6+3 1.332 90 2.929 0.053 12 3.5 0 0.2 3.975 4.005 1.509 -param Pd4+2 1.338 90 2.899 0.048 12 3.21 0 0.2 4.32 4 1.544 -param Ag1+1 1.386 180 3.148 0.036 12 1.956 0 0.2 4.436 3.134 1.622 -param Cd3+2 1.403 109.47 2.848 0.228 12 1.65 0 0.2 5.034 3.957 1.6 -param In3+3 1.459 109.47 4.463 0.599 11 2.07 0 0.2 3.506 2.896 1.404 -param Sn3 1.398 109.47 4.392 0.567 12 2.961 0.199 0.2 3.987 3.124 1.354 -param Sb3+3 1.407 91.6 4.42 0.449 13 2.704 1.1 0.2 4.899 3.342 1.404 -param Te3+2 1.386 90.25 4.47 0.398 14 2.882 0.3 0.2 5.816 3.526 1.38 -param I_ 1.382 180 4.5 0.339 15 2.65 0 0.2 6.822 3.762 1.333 -param Xe4+4 1.267 90 4.404 0.332 12 0.556 0 0.2 7.595 4.975 2.459 -param Cs 2.57 180 4.517 0.045 12 1.573 0 0.1 2.183 1.711 2.984 -param Ba6+2 2.277 90 3.703 0.364 12 2.727 0 0.1 2.814 2.396 2.442 -param La3+3 1.943 109.47 3.522 0.017 12 3.3 0 0.1 2.8355 2.7415 2.071 -param Ce6+3 1.841 90 3.556 0.013 12 3.3 0 0.1 2.774 2.692 1.925 -param Pr6+3 1.823 90 3.606 0.01 12 3.3 0 0.1 2.858 2.564 2.007 -param Nd6+3 1.816 90 3.575 0.01 12 3.3 0 0.1 2.8685 2.6205 2.007 -param Pm6+3 1.801 90 3.547 0.009 12 3.3 0 0.1 2.881 2.673 2 -param Sm6+3 1.78 90 3.52 0.008 12 3.3 0 0.1 2.9115 2.7195 1.978 -param Eu6+3 1.771 90 3.493 0.008 12 3.3 0 0.1 2.8785 2.7875 2.227 -param Gd6+3 1.735 90 3.368 0.009 12 3.3 0 0.1 3.1665 2.9745 1.968 -param Tb6+3 1.732 90 3.451 0.007 12 3.3 0 0.1 3.018 2.834 1.954 -param Dy6+3 1.71 90 3.428 0.007 12 3.3 0 0.1 3.0555 2.8715 1.934 -param Ho6+3 1.696 90 3.409 0.007 12 3.416 0 0.1 3.127 2.891 1.925 -param Er6+3 1.673 90 3.391 0.007 12 3.3 0 0.1 3.1865 2.9145 1.915 -param Tm6+3 1.66 90 3.374 0.006 12 3.3 0 0.1 3.2514 2.9329 2 -param Yb6+3 1.637 90 3.355 0.228 12 2.618 0 0.1 3.2889 2.965 2.158 -param Lu6+3 1.671 90 3.64 0.041 12 3.271 0 0.1 2.9629 2.4629 1.896 -param Hf3+4 1.611 109.47 3.141 0.072 12 3.921 0 0.1 3.7 3.4 1.759 -param Ta3+5 1.511 109.47 3.17 0.081 12 4.075 0 0.1 5.1 2.85 1.605 -param W_6+6 1.392 90 3.069 0.067 12 3.7 0 0.1 4.63 3.31 1.538 -param W_3+4 1.526 109.47 3.069 0.067 12 3.7 0 0.1 4.63 3.31 1.538 -param W_3+6 1.38 109.47 3.069 0.067 12 3.7 0 0.1 4.63 3.31 1.538 -param Re6+5 1.372 90 2.954 0.066 12 3.7 0 0.1 3.96 3.92 1.6 -param Re3+7 1.314 109.47 2.954 0.066 12 3.7 0 0.1 3.96 3.92 1.6 -param Os6+6 1.372 90 3.12 0.037 12 3.7 0 0.1 5.14 3.63 1.7 -param Ir6+3 1.371 90 2.84 0.073 12 3.731 0 0.1 5 4 1.866 -param Pt4+2 1.364 90 2.754 0.08 12 3.382 0 0.1 4.79 4.43 1.557 -param Au4+3 1.262 90 3.293 0.039 12 2.625 0 0.1 4.894 2.586 1.618 -param Hg1+2 1.34 180 2.705 0.385 12 1.75 0 0.1 6.27 4.16 1.6 -param Tl3+3 1.518 120 4.347 0.68 11 2.068 0 0.1 3.2 2.9 1.53 -param Pb3 1.459 109.47 4.297 0.663 12 2.846 0.1 0.1 3.9 3.53 1.444 -param Bi3+3 1.512 90 4.37 0.518 13 2.47 1 0.1 4.69 3.74 1.514 -param Po3+2 1.5 90 4.709 0.325 14 2.33 0.3 0.1 4.21 4.21 1.48 -param At 1.545 180 4.75 0.284 15 2.24 0 0.1 4.75 4.75 1.47 -param Rn4+4 1.42 90 4.765 0.248 16 0.583 0 0.1 5.37 5.37 2.2 -param Fr 2.88 180 4.9 0.05 12 1.847 0 0 2 2 2.3 -param Ra6+2 2.512 90 3.677 0.404 12 2.92 0 0 2.843 2.434 2.2 -param Ac6+3 1.983 90 3.478 0.033 12 3.9 0 0 2.835 2.835 2.108 -param Th6+4 1.721 90 3.396 0.026 12 4.202 0 0 3.175 2.905 2.018 -param Pa6+4 1.711 90 3.424 0.022 12 3.9 0 0 2.985 2.905 1.8 -param U_6+4 1.684 90 3.395 0.022 12 3.9 0 0 3.341 2.853 1.713 -param Np6+4 1.666 90 3.424 0.019 12 3.9 0 0 3.549 2.717 1.8 -param Pu6+4 1.657 90 3.424 0.016 12 3.9 0 0 3.243 2.819 1.84 -param Am6+4 1.66 90 3.381 0.014 12 3.9 0 0 2.9895 3.0035 1.942 -param Cm6+3 1.801 90 3.326 0.013 12 3.9 0 0 2.8315 3.1895 1.9 -param Bk6+3 1.761 90 3.339 0.013 12 3.9 0 0 3.1935 3.0355 1.9 -param Cf6+3 1.75 90 3.313 0.013 12 3.9 0 0 3.197 3.101 1.9 -param Es6+3 1.724 90 3.299 0.012 12 3.9 0 0 3.333 3.089 1.9 -param Fm6+3 1.712 90 3.286 0.012 12 3.9 0 0 3.4 3.1 1.9 -param Md6+3 1.689 90 3.274 0.011 12 3.9 0 0 3.47 3.11 1.9 -param No6+3 1.679 90 3.248 0.011 12 3.9 0 0 3.475 3.175 1.9 -param Lw6+3 1.698 90 3.236 0.011 12 3.9 0 0 3.5 3.2 1.9 This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2012-05-14 03:17:26
|
Revision: 17126 http://jmol.svn.sourceforge.net/jmol/?rev=17126&view=rev Author: hansonr Date: 2012-05-14 03:17:17 +0000 (Mon, 14 May 2012) Log Message: ----------- MMFF94 done: bond-stretch, angle-bend, torsion-twist. out-of-plane bending. Still to do: electrostatic charges, stretch-bend, and van der Waals. Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/MinAtom.java trunk/Jmol/src/org/jmol/minimize/Minimizer.java trunk/Jmol/src/org/jmol/minimize/Util.java trunk/Jmol/src/org/jmol/minimize/forcefield/Calculation.java trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java Modified: trunk/Jmol/src/org/jmol/minimize/MinAtom.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/MinAtom.java 2012-05-13 04:29:32 UTC (rev 17125) +++ trunk/Jmol/src/org/jmol/minimize/MinAtom.java 2012-05-14 03:17:17 UTC (rev 17126) @@ -32,30 +32,31 @@ public class MinAtom { - int index; + int rawIndex; public String sType; public Atom atom; public AtomType ffAtomType; public int ffType; + public Integer vdwKey; public double[] coord = new double[3]; public double[] force = new double[3]; public List<MinBond> bonds = new ArrayList<MinBond>(); public int nBonds; public int hCount; public double partialCharge; - + int[] bondedAtoms; - + @Override public String toString() { - return (sType == null ? atom.getAtomName() : "#" + index + " " + sType); + return "#" + rawIndex + " " + sType; } - + MinAtom(int index, Atom atom, double[] coord) { - this.index = index; + this.rawIndex = index; this.atom = atom; this.coord = coord; - hCount = atom.getCovalentHydrogenCount(); + hCount = atom.getCovalentHydrogenCount(); } void set() { @@ -76,7 +77,7 @@ if (bondedAtoms == null) { bondedAtoms = new int[nBonds]; for (int i = nBonds; --i >= 0;) - bondedAtoms[i] = bonds.get(i).getOtherAtom(index); + bondedAtoms[i] = bonds.get(i).getOtherAtom(rawIndex); } return bondedAtoms; } Modified: trunk/Jmol/src/org/jmol/minimize/Minimizer.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2012-05-13 04:29:32 UTC (rev 17125) +++ trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2012-05-14 03:17:17 UTC (rev 17126) @@ -288,6 +288,7 @@ bsElements.set(atomicNo); minAtoms[pt] = new MinAtom(pt, atom, new double[] { atom.x, atom.y, atom.z }); + minAtoms[pt].sType = atom.getAtomName(); } Logger.info(GT._("{0} atoms will be minimized.", "" + atomCount)); @@ -306,12 +307,24 @@ return true; } + private void setAtomPositions() { + for (int i = 0; i < atomCount; i++) + minAtoms[i].set(); + bsMinFixed = null; + if (bsFixed != null) { + bsMinFixed = new BitSet(); + for (int i = bsAtoms.nextSetBit(0), pt = 0; i >= 0; i = bsAtoms + .nextSetBit(i + 1), pt++) + if (bsFixed.get(i)) + bsMinFixed.set(pt); + } + } + private void getBonds() { - // add all bonds List<MinBond> bondInfo = new ArrayList<MinBond>(); bondCount = 0; int i1, i2; - for (int i = rawBondCount; --i >= 0;) { + for (int i = 0; i < rawBondCount; i++) { Bond bond = bonds[i]; if (!bsAtoms.get(i1 = bond.getAtomIndex1()) || !bsAtoms.get(i2 = bond.getAtomIndex2())) @@ -349,21 +362,7 @@ minAtoms[i].getBondedAtomIndexes(); } - private void setAtomPositions() { - for (int i = 0; i < atomCount; i++) - minAtoms[i].set(); - bsMinFixed = null; - if (bsFixed != null) { - bsMinFixed = new BitSet(); - for (int i = bsAtoms.nextSetBit(0), pt = 0; i >= 0; i = bsAtoms - .nextSetBit(i + 1), pt++) - if (bsFixed.get(i)) - bsMinFixed.set(pt); - } - } - public void getAngles() { - List<MinAngle> vAngles = new ArrayList<MinAngle>(); int[] atomList; int ic; @@ -406,7 +405,7 @@ List<MinTorsion> vTorsions = new ArrayList<MinTorsion>(); int id; // extend all angles a-b-c by one, but only - // when when c > b or b > a + // when when c > b or a > b for (int i = minAngles.length; --i >= 0;) { int[] angle = minAngles[i].data; int ia = angle[0]; @@ -449,10 +448,7 @@ minTorsions = vTorsions.toArray(new MinTorsion[vTorsions.size()]); Logger.info(minTorsions.length + " torsions"); } - - - ///////////////////////////// minimize ////////////////////// public ForceField getForceField(String ff) { Modified: trunk/Jmol/src/org/jmol/minimize/Util.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Util.java 2012-05-13 04:29:32 UTC (rev 17125) +++ trunk/Jmol/src/org/jmol/minimize/Util.java 2012-05-14 03:17:17 UTC (rev 17126) @@ -116,23 +116,41 @@ return (canBeSquared(a.x) && canBeSquared(a.y) && canBeSquared(a.z)); } - /* Calculate the angle between point a and the plane determined by b,c,d */ + /** + * + * calculates angle of a to plane bcd, returning a value > pi/2 in + * highly distorted trigonal pyramidal situations + * + * @param a + * @param b + * @param c + * @param d + * @param v1 + * @param v2 + * @param norm + * @return Wilson angle + */ public static double pointPlaneAngleRadians(Vector3d a, Vector3d b, Vector3d c, Vector3d d, Vector3d v1,Vector3d v2, Vector3d norm) { + v1.sub(b, c); v2.sub(b, d); - norm.cross(v1, v2); + norm.cross(v1, v2); v2.add(v1); v1.sub(b, a); - double angleA_CD = vectorAngleRadians(v2, v1); + // we need to allow for very distorted cases, where a, c, and d are very close to each other + double angleA_CD = vectorAngleRadians(v2, v1); + // normally angleA_CD is > pi/2, because v1 is from a to b double angleNorm = vectorAngleRadians(norm, v1); + // angleNorm could be > PI/2, so we correct for that here: if (angleNorm > Math.PI / 2) angleNorm = Math.PI - angleNorm; - return Math.PI / 2.0 + (angleA_CD > Math.PI / 2.0 ? -angleNorm : angleNorm) ; + // for highly distorted groups, return will be > pi/2 + return Math.PI / 2.0 + (angleA_CD > Math.PI / 2.0 ? -angleNorm : angleNorm) ; } - + private static double vectorAngleRadians(Vector3d v1, Vector3d v2) { double l1 = v1.length(); double l2 = v2.length(); Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/Calculation.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/Calculation.java 2012-05-13 04:29:32 UTC (rev 17125) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/Calculation.java 2012-05-14 03:17:17 UTC (rev 17126) @@ -46,7 +46,18 @@ } void getPointers(Object[] dataIn) { + dData = (double[])dataIn[1]; iData = (int[])dataIn[0]; - dData = (double[])dataIn[1]; + switch (iData.length) { + default: + id = iData[3]; + case 3: + ic = iData[2]; + case 2: + ib = iData[1]; + case 1: + ia = iData[0]; + case 0: + } } } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java 2012-05-13 04:29:32 UTC (rev 17125) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java 2012-05-14 03:17:17 UTC (rev 17126) @@ -57,6 +57,8 @@ int atomCount; int bondCount; + int angleCount; + int torsionCount; MinAtom[] minAtoms; MinBond[] minBonds; MinAngle[] minAngles; @@ -80,9 +82,11 @@ this.minBonds = minBonds; this.minAngles = minAngles; this.minTorsions = minTorsions; - this.constraints = constraints; atomCount = minAtoms.length; bondCount = minBonds.length; + angleCount = minAngles.length; + torsionCount = minTorsions.length; + this.constraints = constraints; if (partialCharges != null && partialCharges.length == atomCount) for (int i = atomCount; --i >= 0;) if (partialCharges[i] != 0) { @@ -96,12 +100,6 @@ abstract boolean setupCalculations(); - abstract String getAtomList(String title); - - abstract String getDebugHeader(int iType); - - abstract String getDebugFooter(int iType, double energy); - abstract String getUnit(); abstract double compute(int iType, Object[] dataIn); @@ -127,7 +125,6 @@ } void appendLogData(String s) { - System.out.println("LOG " + s); logData.append(s).append("\n"); } @@ -404,4 +401,246 @@ appendLogData("---------------------\n"); } + String getAtomList(String title) { + String trailer = + "----------------------------------------" + + "-------------------------------------------------------\n"; + StringBuffer sb = new StringBuffer(); + sb.append("\n" + title + "\n\n" + + " ATOM X Y Z TYPE GRADX GRADY GRADZ " + + "---------BONDED ATOMS--------\n" + + trailer); + for (int i = 0; i < atomCount; i++) { + MinAtom atom = minAtoms[i]; + int[] others = atom.getBondedAtomIndexes(); + int[] iVal = new int[others.length + 1]; + iVal[0] = atom.atom.getAtomNumber(); + String s = " "; + for (int j = 0; j < others.length; j++) { + s += " %3d"; + iVal[j + 1] = minAtoms[others[j]].atom.getAtomNumber(); + } + sb.append(TextFormat.sprintf("%3d %8.3f %8.3f %8.3f %-5s %8.3f %8.3f %8.3f" + s + "\n", + new Object[] { atom.sType, + new float[] { (float) atom.coord[0], (float) atom.coord[1], + (float) atom.coord[2], (float) atom.force[0], (float) atom.force[1], + (float) atom.force[2] }, iVal})); + } + sb.append(trailer + "\n\n"); + return sb.toString(); + } + + String getDebugHeader(int iType) { + switch (iType){ + case -1: + //Override to give reference + break; + case CALC_DISTANCE: + return + "\nB O N D S T R E T C H I N G (" + bondCount + " bonds)\n\n" + +" ATOMS ATOM TYPES BOND BOND IDEAL FORCE\n" + +" I J I J TYPE LENGTH LENGTH CONSTANT DELTA ENERGY\n" + +"--------------------------------------------------------------------------------"; + case CALC_ANGLE: + return + "\nA N G L E B E N D I N G (" + minAngles.length + " angles)\n\n" + +" ATOMS ATOM TYPES VALENCE IDEAL FORCE\n" + +" I J K I J K ANGLE ANGLE CONSTANT ENERGY\n" + +"--------------------------------------------------------------------------"; + case CALC_TORSION: + return + "\nT O R S I O N A L (" + minTorsions.length + " torsions)\n\n" + +" ATOMS ATOM TYPES n COS FORCE TORSION\n" + +" I J K L I J K L (n phi0) CONSTANT ANGLE ENERGY\n" + +"---------------------------------------------------------------------------------------------"; + case CALC_OOP: + return + "\nO U T - O F - P L A N E B E N D I N G\n\n" + +" ATOMS ATOM TYPES OOP FORCE \n" + +" I J K L I J K L ANGLE CONSTANT ENERGY\n" + +"--------------------------------------------------------------------------"; + case CALC_VDW: + return + "\nV A N D E R W A A L S\n\n" + +" ATOMS ATOM TYPES\n" + +" I J I J Rij kij ENERGY\n" + +"-----------------------------------------------"; + case CALC_ES: + return + "\nE L E C T R O S T A T I C I N T E R A C T I O N S\n\n" + +" ATOMS ATOM TYPES QiQj\n" + +" I J I J Rij *332.17 ENERGY\n" + +"-----------------------------------------------"; + } + return ""; + } + + String getDebugLine(int iType, Calculation c) { + switch (iType) { + case CALC_DISTANCE: + return TextFormat.sprintf( + "%3d %3d %-5s %-5s %4.2f%8.3f %8.3f %8.3f %8.3f %8.3f", + new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, + new float[] { 0, (float)c.rab, + (float)c.dData[1], (float)c.dData[0], + (float)c.delta, (float)c.energy }, + new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() }}); + case CALC_ANGLE: + return TextFormat.sprintf( + "%3d %3d %3d %-5s %-5s %-5s %8.3f %8.3f %8.3f %8.3f", + new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, + minAtoms[c.ic].sType, + new float[] { (float)(c.theta * RAD_TO_DEG), (float) c.dData[1] /*THETA0*/, + (float)c.dData[0]/*Kijk*/, (float) c.energy }, + new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), + minAtoms[c.ic].atom.getAtomNumber()} }); + case CALC_TORSION: + return TextFormat.sprintf( + "%3d %3d %3d %3d %-5s %-5s %-5s %-5s %3d %8.3f %8.3f %8.3f %8.3f", + new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, + minAtoms[c.ic].sType, minAtoms[c.id].sType, + new float[] { (float) c.dData[1]/*cosNphi0*/, (float) c.dData[0]/*V*/, + (float) (c.theta * RAD_TO_DEG), (float) c.energy }, + new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), + minAtoms[c.ic].atom.getAtomNumber(), minAtoms[c.id].atom.getAtomNumber(), c.iData[4] } }); + case CALC_OOP: + return TextFormat.sprintf("" + + "%3d %3d %3d %3d %-5s %-5s %-5s %-5s %8.3f %8.3f %8.3f", + new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, + minAtoms[c.ic].sType, minAtoms[c.id].sType, + new float[] { (float)(c.theta * RAD_TO_DEG), + (float)c.dData[0]/*koop*/, (float) c.energy }, + new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), + minAtoms[c.ic].atom.getAtomNumber(), minAtoms[c.id].atom.getAtomNumber() } }); + case CALC_VDW: + return TextFormat.sprintf("%3d %3d %-5s %-5s %6.3f %8.3f %8.3f", + new Object[] { minAtoms[c.iData[0]].sType, minAtoms[c.iData[1]].sType, + new float[] { (float)c.rab, (float)c.dData[0]/*kab*/, (float)c.energy}, + new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() } }); + case CALC_ES: + return TextFormat.sprintf("%3d %3d %-5s %-5s %6.3f %8.3f %8.3f", + new Object[] { minAtoms[c.iData[0]].sType, minAtoms[c.iData[1]].sType, + new float[] { (float)c.rab, (float)c.dData[0]/*qq*/, (float)c.energy }, + new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() } }); + } + return ""; + } + + String getDebugFooter(int iType, double energy) { + String s = ""; + switch (iType){ + case CALC_DISTANCE: + s = "BOND STRETCHING"; + break; + case CALC_ANGLE: + s = "ANGLE BENDING"; + break; + case CALC_TORSION: + s = "TORSIONAL"; + break; + case CALC_OOP: + s = "OUT-OF-PLANE BENDING"; + break; + case CALC_VDW: + s = "VAN DER WAALS"; + break; + case CALC_ES: + s = "ELECTROSTATIC ENERGY"; + break; + } + return TextFormat.sprintf("\n TOTAL %s ENERGY = %8.3f %s\n", + new Object[] { s, getUnit(), Float.valueOf((float) energy) }); + } + + //////////////////////////////////////////////////// + + void setBondVariables(Calculation c) { + if (gradients) { + setCoords(c, 2); + c.rab = Util.restorativeForceAndDistance(da, db, dc); + } else { + c.rab = Math.sqrt(Util.distance2(minAtoms[c.ia].coord, minAtoms[c.ib].coord)); + } + } + + void setAngleVariables(Calculation c) { + if (gradients) { + setCoords(c, 3); + c.theta = Util.restorativeForceAndAngleRadians(da, db, dc); + } else { + c.theta = Util.getAngleRadiansABC(minAtoms[c.ia].coord, minAtoms[c.ib].coord, minAtoms[c.ic].coord); + } + if (!Util.isFinite(c.theta)) + c.theta = 0.0; + } + + void setOopVariables(Calculation c) { + setCoords(c, 4); + if (gradients) { + c.theta = Util.restorativeForceAndOutOfPlaneAngleRadians(da, db, dc, dd, v1, v2, v3); + } else { + //if (atoms[ib].atom.getAtomIndex() == 20) + //System.out.println(" atom palne angl " + atoms[ib].atom.getInfo()); + + c.theta = Util.pointPlaneAngleRadians(da, db, dc, dd, v1, v2, v3); + } + if (!Util.isFinite(c.theta)) + c.theta = 0.0; + } + + void setTorsionVariables(Calculation c) { + if (gradients) { + setCoords(c, 4); + c.theta = Util.restorativeForceAndTorsionAngleRadians(da, db, dc, dd); + if (!Util.isFinite(c.theta)) + c.theta = 0.001 * DEG_TO_RAD; + } else { + c.theta = Util.getTorsionAngleRadians(minAtoms[c.ia].coord, minAtoms[c.ib].coord, + minAtoms[c.ic].coord, minAtoms[c.id].coord, v1, v2, v3); + } + } + + void setPairVariables(Calculation c) { + if (gradients) { + setCoords(c, 2); + c.rab = Util.restorativeForceAndDistance(da, db, dc); + } else { + c.rab = Math.sqrt(Util.distance2(minAtoms[c.ia].coord, minAtoms[c.ib].coord)); + } + if (Util.isNearZero(c.rab, 1.0e-3)) + c.rab = 1.0e-3; + } + + private void setCoords(Calculation c, int n) { + switch(n) { + case 4: + da.set(minAtoms[c.ia].coord); + // fall through + case 3: + db.set(minAtoms[c.ib].coord); + // fall through + case 2: + dc.set(minAtoms[c.ic].coord); + // fall through + case 1: + dd.set(minAtoms[c.id].coord); + } + } + + void addForces(Calculation c, int n) { + switch (n) { + case 4: + addForce(dd, c.id, c.dE); + // fall through + case 3: + addForce(dc, c.ic, c.dE); + // fall through + case 2: + addForce(db, c.ib, c.dE); + // fall through + case 1: + addForce(da, c.ia, c.dE); + } + } + } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java 2012-05-13 04:29:32 UTC (rev 17125) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java 2012-05-14 03:17:17 UTC (rev 17126) @@ -65,8 +65,11 @@ class CalculationsMMFF extends Calculations { - public static final double fStretch = 143.9325 / 2; + final static double FPAR = 143.9325; + public static final int DA_D = 'D'; + public static final int DA_DA = DA_D + 'A'; + private Map<Integer, Object> ffParams; DistanceCalc bondCalc; @@ -109,28 +112,25 @@ calc = calculations[CALC_ANGLE] = new ArrayList<Object[]>(); AngleCalc angleCalc = new AngleCalc(); - for (int i = minAngles.length; --i >= 0;) - angleCalc.setData(calc, minAngles[i].data); + for (int i = 0; i < angleCount; i++) + angleCalc.setData(calc, minAngles[i]); calc = calculations[CALC_STRETCH_BEND] = new ArrayList<Object[]>(); SBCalc sbCalc = new SBCalc(); - for (int i = minAngles.length; --i >= 0;) - sbCalc.setData(calc, minAngles[i].data); + for (int i = 0; i < angleCount; i++) + sbCalc.setData(calc, minAngles[i]); calc = calculations[CALC_TORSION] = new ArrayList<Object[]>(); TorsionCalc torsionCalc = new TorsionCalc(); - for (int i = minTorsions.length; --i >= 0;) - torsionCalc.setData(calc, minTorsions[i].data); + for (int i = 0; i < torsionCount; i++) + torsionCalc.setData(calc, minTorsions[i]); calc = calculations[CALC_OOP] = new ArrayList<Object[]>(); // set up the special atom arrays OOPCalc oopCalc = new OOPCalc(); - int elemNo; - for (int i = 0; i < atomCount; i++) { - MinAtom a = minAtoms[i]; - if (a.nBonds == 3 && isInvertible(elemNo = a.atom.getElementNumber())) - oopCalc.setData(calc, i, elemNo); - } + for (int i = 0; i < atomCount; i++) + if (isInvertible(minAtoms[i])) + oopCalc.setData(calc, i); pairSearch(calculations[CALC_VDW] = new ArrayList<Object[]>(), new VDWCalc(), calculations[CALC_ES] = new ArrayList<Object[]>(), new ESCalc()); @@ -138,36 +138,45 @@ return true; } - private boolean isInvertible(int n) { - switch (n) { - case 6: // C - case 7: // N - case 8: // O - case 15: // P - case 33: // As - case 51: // Sb - case 83: // Bi + private static boolean isInvertible(MinAtom a) { + + // defined for typeB = 2, 3, 8, 10, 17, 26, 30, 37, 39, 40, 41, 43, + // 45, 49, 54, 55, 56, 57, 58, 63, 64, 67, 69, 73, 78, 80, 81, 82 + // but is 0 for 8 (amines), 17 (sulfoxides), + // 26 (PD3), 43 (N-S), 73 (O-S(=O)R, 82 (N-oxide) + // that is, just the planar ones: + // 2, 3, 10, 30, 37, 39, 40, 41, + // 45, 49, 54, 55, 56, 57, 58, 63, + // 64, 67, 69, 78, 80, 81 + switch (a.ffType) { + default: + return false; + case 2: + case 3: + case 10: + case 30: + case 37: + case 39: + case 40: + case 41: + case 45: + case 49: + case 54: + case 55: + case 56: + case 57: + case 58: + case 63: + case 64: + case 67: + case 69: + case 78: + case 80: + case 81: return true; - default: - return false;// no inversion term for this element } } - static double calculateR0(double ri, double rj, double chiI, double chiJ, - double bondorder) { - // precompute the equilibrium geometry - // From equation 3 - double rbo = -0.1332 * (ri + rj) * Math.log(bondorder); - // From equation 4 - - double dchi = Math.sqrt(chiI) - Math.sqrt(chiJ); - double ren = ri * rj * dchi * dchi / (chiI * ri + chiJ * rj); - // From equation 2 - // NOTE: See http://towhee.sourceforge.net/forcefields/uff.html - // There is a typo in the published paper - return (ri + rj + rbo - ren); - } - @Override double compute(int iType, Object[] dataIn) { @@ -194,52 +203,42 @@ return ffParams.get(a); } - double r0, kb; - final static double cs = -2.0; - final static double cs2 = ((7.0/12.0)*(cs * cs)); - double delta2; class DistanceCalc extends Calculation { + final static double FSTRETCH = FPAR / 2; + final static double CS = -2.0; + final static double CS2 = ((7.0/12.0)*(CS * CS)); + + double r0, kb; + double delta2; + void setData(List<Object[]> calc, MinBond bond) { ia = bond.data[0]; ib = bond.data[1]; - - double[] data = (double[]) getParameter(bond.key); -/* - - kb = KCAL332 * parA.dVal[PAR_Z] * parB.dVal[PAR_Z] / (r0 * r0 * r0); -*/ + Object data = getParameter(bond.key); + if (data == null) + return; calc.add(new Object[] { new int[] { ia, ib }, data }); } @Override double compute(Object[] dataIn) { + getPointers(dataIn); kb = dData[0]; r0 = dData[1]; - ia = iData[0]; - ib = iData[1]; + setBondVariables(this); - if (gradients) { - da.set(minAtoms[ia].coord); - db.set(minAtoms[ib].coord); - rab = Util.restorativeForceAndDistance(da, db, dc); - } else { - rab = Math.sqrt(Util.distance2(minAtoms[ia].coord, minAtoms[ib].coord)); - } - - // Er = 0.5 k (r - r0)^2 - - delta = rab - r0; // we pre-compute the r0 below + delta = rab - r0; delta2 = delta * delta; - energy = kb * delta * delta; // 0.5 factor was precalculated - energy = fStretch * kb * (delta2) * (1 + cs * delta + cs2 * (delta2)); + energy = FSTRETCH * kb * delta2 * (1 + CS * delta + CS2 * (delta2)); if (gradients) { - dE = 2.0 * kb * delta; // TODO - addForce(da, ia, dE); - addForce(db, ib, dE); + //TODO CHECK THIS CAREFULLY + dE = FSTRETCH * 2 * kb * delta * (1 + CS * delta + CS2 * delta2 + + 0.5 * delta * (CS + 2 * CS2 * delta)); + addForces(this, 2); } if (logging) @@ -250,23 +249,53 @@ } - final static double KCAL644 = 644.12 * KCAL_TO_KJ; - class AngleCalc extends Calculation { - - void setData(List<Object[]> calc, int[] angle) { + + void setData(List<Object[]> calc, MinAngle angle) { + Object data = getParameter(angle.key); + if (data == null) + return; + calc.add(new Object[] { angle.data, data }); } + final static double CB = -0.4 * DEG_TO_RAD; + @Override double compute(Object[] dataIn) { + + getPointers(dataIn); + double ka = dData[0]; + double t0 = dData[1]; + setAngleVariables(this); + + double dt = (theta * RAD_TO_DEG - t0); + + // could have problems here for very distorted structures. + + if (t0 == 180) { + energy = FPAR * ka * (1 + Math.cos(theta)); + if (gradients) + dE = -FPAR * ka * Math.sin(theta); + } else { + energy = 0.021922 * ka * Math.pow(dt, 2) * (1 + CB * dt); // 0.043844/2 + if (gradients) + dE = 0.021922 * ka * dt * (2 + 3 * CB * dt); + } + if (gradients) + addForces(this, 3); + + if (logging) + appendLogData(getDebugLine(CALC_ANGLE, this)); + return energy; } + } class SBCalc extends Calculation { // TODO - void setData(List<Object[]> calc, int[] angle) { + void setData(List<Object[]> calc, MinAngle angle) { } @Override @@ -277,26 +306,98 @@ class TorsionCalc extends Calculation { - void setData(List<Object[]> calc, int[] t) { + void setData(List<Object[]> calc, MinTorsion t) { + Object data = getParameter(t.key); + if (data == null) + return; + calc.add(new Object[] { t.data, data }); } @Override double compute(Object[] dataIn) { + + getPointers(dataIn); + double v1 = dData[0]; + double v2 = dData[1]; + double v3 = dData[2]; + setTorsionVariables(this); + + // use one single cosine calculation + + double cosTheta = Math.cos(theta); + double cosTheta2 = cosTheta * cosTheta; + + energy = 0.5 * (v1 * (1 + cosTheta) + + v2 * (2 - 2 * cosTheta2) + + v3 * (1 + cosTheta * (4 * cosTheta2 - 3))); + +/* + energy = 0.5 * (v1 * (1.0 + Math.cos(theta)) + + v2 * (1 - Math.cos(2 * theta)) + + v3 * (1 + Math.cos(3 * theta))); +*/ + if (gradients) { + double sinTheta = Math.sin(theta); + dE = 0.5 * (-v1 * sinTheta + + 4 * v2 * sinTheta * cosTheta + + 3 * v3 * sinTheta * (1 - 4 * cosTheta2)); +/* + dE = 0.5 * (-v1 * sinTheta + + 2 * v2 * Math.sin(2 * theta) + - 3 * v3 * Math.sin(3 * theta)); +*/ + addForces(this, 4); + } + + if (logging) + appendLogData(getDebugLine(CALC_TORSION, this)); + return energy; } } - final static double KCAL6 = 6.0 * KCAL_TO_KJ; - final static double KCAL22 = 22.0 * KCAL_TO_KJ; - final static double KCAL44 = 44.0 * KCAL_TO_KJ; - class OOPCalc extends Calculation { - void setData(List<Object[]> calc, int ib, int elemNo) { + final static double FOOPD = 0.043844 * RAD_TO_DEG; + final static double FOOP = FOOPD / 2 * RAD_TO_DEG; + + int[] list = new int[4]; + + void setData(List<Object[]> calc, int i) { + if (minAtoms[i].nBonds != 3) + return;// should not be possible... + int[] indices = minAtoms[i].getBondedAtomIndexes(); + // our calculation is for first, not last, relative to plane of others, + list[0] = indices[2]; + list[1] = i; + list[2] = indices[1]; + list[3] = indices[0]; + double koop = ((ForceFieldMMFF) ff).getOutOfPlaneParameter(list); + if (koop == 0) + return; + double[] dk = new double[] { koop }; + calc.add(new Object[] { new int[] { indices[0], i, indices[1], indices[2] }, dk }); + calc.add(new Object[] { new int[] { indices[1], i, indices[2], indices[0] }, dk }); + calc.add(new Object[] { new int[] { indices[2], i, indices[0], indices[1] }, dk }); } @Override double compute(Object[] dataIn) { + + getPointers(dataIn); + setOopVariables(this); + double koop = dData[0]; + + energy = FOOP * koop * theta * theta; // theta in radians + + if (gradients) { + dE = FOOPD * koop * theta; + addForces(this, 4); + } + + if (logging) + appendLogData(getDebugLine(CALC_OOP, this)); + return energy; } } @@ -305,16 +406,67 @@ @Override void setData(List<Object[]> calc, int ia, int ib) { + a = minAtoms[ia]; + b = minAtoms[ib]; + double[] dataA = (double[]) getParameter(a.vdwKey); + double[] dataB = (double[]) getParameter(a.vdwKey); + if (dataA == null || dataB == null) + return; + double alpha_a = dataA[0]; + double N_a = dataA[1]; + double A_a = dataA[2]; + double G_a = dataA[3]; + int DA_a = (int) dataA[4]; + double alpha_b = dataB[0]; + double N_b = dataB[1]; + double A_b = dataB[2]; + double G_b = dataB[3]; + int DA_b = (int) dataB[4]; + + + double rs_aa = A_a * Math.pow(alpha_a, 0.25); + double rs_bb = A_b * Math.pow(alpha_b, 0.25); + double gamma = (rs_aa - rs_bb) / (rs_aa + rs_bb); + double rs = 0.5 * (rs_aa + rs_bb); + if (DA_a != DA_D && DA_b != DA_D) + rs *= (1.0 + 0.2 * (1.0 - Math.exp(-12.0 * gamma * gamma))); + double eps = ((181.16 * G_a * G_b * alpha_a * alpha_b) + / (Math.sqrt(alpha_a / N_a) + Math.sqrt(alpha_b / N_b))) * Math.pow(rs, -6.0); + + if(DA_a + DA_b == DA_DA) { + rs *= 0.8; + eps *= 0.5; + } + calc.add(new Object[] { new int[] {ia, ib}, new double[] { rs, eps } }); } @Override double compute(Object[] dataIn) { + getPointers(dataIn); + double rs = dData[0]; + double eps = dData[1]; + + rab = Math.sqrt(Util.distance2(minAtoms[ia].coord, minAtoms[ib].coord)); + double r_rs = rab / rs; + energy = eps * Math.pow(1.07/ (r_rs + 0.07), 7) * (1.12 / (Math.pow(r_rs, 7) + 0.12) - 2); + + if (gradients) { + //TODO CHECK CAREFULLY + dE = 7 * eps * Math.pow(1.07 / (r_rs + 0.07), 6) * + ((-1.07 * rs / Math.pow(rab + 0.07 * rs, 2)) + * (1.12 / (Math.pow(r_rs, 7) + 0.12) - 2) + + (-1.12 * Math.pow(rs, 7) / rab + + 0.12 * Math.pow(rs, 14)) * (1.07 / (r_rs + 0.07))); + addForces(this, 2); + } + + if (logging) + appendLogData(getDebugLine(CALC_VDW, this)); + return energy; } } - final static double KCAL332 = KCAL_TO_KJ * 332.0637; - class ESCalc extends PairCalc { @Override @@ -326,163 +478,43 @@ return energy; } } - ///////// REPORTING ///////////// @Override - String getAtomList(String title) { - String trailer = - "----------------------------------------" - + "-------------------------------------------------------\n"; - StringBuffer sb = new StringBuffer(); - sb.append("\n" + title + "\n\n" - + " ATOM X Y Z TYPE GRADX GRADY GRADZ " - + "---------BONDED ATOMS--------\n" - + trailer); - for (int i = 0; i < atomCount; i++) { - MinAtom atom = minAtoms[i]; - int[] others = atom.getBondedAtomIndexes(); - int[] iVal = new int[others.length + 1]; - iVal[0] = atom.atom.getAtomNumber(); - String s = " "; - for (int j = 0; j < others.length; j++) { - s += " %3d"; - iVal[j + 1] = minAtoms[others[j]].atom.getAtomNumber(); - } - sb.append(TextFormat.sprintf("%3d %8.3f %8.3f %8.3f %-5s %8.3f %8.3f %8.3f" + s + "\n", - new Object[] { atom.sType, - new float[] { (float) atom.coord[0], (float) atom.coord[1], - (float) atom.coord[2], (float) atom.force[0], (float) atom.force[1], - (float) atom.force[2] }, iVal})); - } - sb.append(trailer + "\n\n"); - return sb.toString(); - } - - @Override String getDebugHeader(int iType) { switch (iType){ case -1: - return "Universal Force Field -- " + - "Rappe, A. K., et. al.; J. Am. Chem. Soc. (1992) 114(25) p. 10024-10035\n"; - case CALC_DISTANCE: - return - "\nB O N D S T R E T C H I N G (" + bondCount + " bonds)\n\n" - +" ATOMS ATOM TYPES BOND BOND IDEAL FORCE\n" - +" I J I J TYPE LENGTH LENGTH CONSTANT DELTA ENERGY\n" - +"--------------------------------------------------------------------------------"; - case CALC_ANGLE: - return - "\nA N G L E B E N D I N G (" + minAngles.length + " angles)\n\n" - +" ATOMS ATOM TYPES VALENCE IDEAL FORCE\n" - +" I J K I J K ANGLE ANGLE CONSTANT ENERGY\n" - +"--------------------------------------------------------------------------"; + return "MMFF94 Force Field -- " + + "T. A. Halgren, J. Comp. Chem. 5 & 6 490-519ff (1996).\n"; case CALC_TORSION: return "\nT O R S I O N A L (" + minTorsions.length + " torsions)\n\n" - +" ATOMS ATOM TYPES n COS FORCE TORSION\n" - +" I J K L I J K L (n phi0) CONSTANT ANGLE ENERGY\n" - +"---------------------------------------------------------------------------------------------"; - case CALC_OOP: - return - "\nO U T - O F - P L A N E B E N D I N G\n\n" - +" ATOMS ATOM TYPES OOP FORCE \n" - +" I J K L I J K L ANGLE CONSTANT ENERGY\n" - +"--------------------------------------------------------------------------"; - case CALC_VDW: - return - "\nV A N D E R W A A L S\n\n" - +" ATOMS ATOM TYPES\n" - +" I J I J Rij kij ENERGY\n" - +"-----------------------------------------------"; - case CALC_ES: - return - "\nE L E C T R O S T A T I C I N T E R A C T I O N S\n\n" - +" ATOMS ATOM TYPES QiQj\n" - +" I J I J Rij *332.17 ENERGY\n" - +"-----------------------------------------------"; + +" ATOMS ATOM TYPES TORSION\n" + +" I J K L I J K L ANGLE V1 V2 V3 ENERGY\n" + +"--------------------------------------------------------------------------------------\n"; + default: + return super.getDebugHeader(iType); } - return ""; } + @Override String getDebugLine(int iType, Calculation c) { switch (iType) { - case CALC_DISTANCE: - return TextFormat.sprintf( - "%3d %3d %-5s %-5s %4.2f%8.3f %8.3f %8.3f %8.3f %8.3f", - new Object[] { minAtoms[c.ia].toString(), minAtoms[c.ib].toString(), - new float[] { 0, (float)c.rab, - (float)c.dData[1], (float)c.dData[0], - (float)c.delta, (float)c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() }}); - case CALC_ANGLE: - return TextFormat.sprintf( - "%3d %3d %3d %-5s %-5s %-5s %8.3f %8.3f %8.3f %8.3f", - new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, - minAtoms[c.ic].sType, - new float[] { (float)(c.theta * RAD_TO_DEG), (float)c.dData[4] /*THETA0*/, - (float)c.dData[0]/*Kijk*/, (float) c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), - minAtoms[c.ic].atom.getAtomNumber()} }); case CALC_TORSION: - return TextFormat.sprintf( - "%3d %3d %3d %3d %-5s %-5s %-5s %-5s %3d %8.3f %8.3f %8.3f %8.3f", - new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, + return TextFormat.sprintf( + "%3d %3d %3d %3d %-5s %-5s %-5s %-5s %8.3f %8.3f %8.3f %8.3f %8.3f", + new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, minAtoms[c.ic].sType, minAtoms[c.id].sType, - new float[] { (float) c.dData[1]/*cosNphi0*/, (float) c.dData[0]/*V*/, - (float) (c.theta * RAD_TO_DEG), (float) c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), - minAtoms[c.ic].atom.getAtomNumber(), minAtoms[c.id].atom.getAtomNumber(), c.iData[4] } }); - case CALC_OOP: - return TextFormat.sprintf("" + - "%3d %3d %3d %3d %-5s %-5s %-5s %-5s %8.3f %8.3f %8.3f", - new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, - minAtoms[c.ic].sType, minAtoms[c.id].sType, - new float[] { (float)(c.theta * RAD_TO_DEG), - (float)c.dData[0]/*koop*/, (float) c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), + new float[] { (float) (c.theta * RAD_TO_DEG), (float) c.dData[0]/*v1*/, (float) c.dData[1]/*v2*/, (float) c.dData[2]/*v3*/, + (float) c.energy }, + new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), minAtoms[c.ic].atom.getAtomNumber(), minAtoms[c.id].atom.getAtomNumber() } }); - case CALC_VDW: - return TextFormat.sprintf("%3d %3d %-5s %-5s %6.3f %8.3f %8.3f", - new Object[] { minAtoms[c.iData[0]].sType, minAtoms[c.iData[1]].sType, - new float[] { (float)c.rab, (float)c.dData[0]/*kab*/, (float)c.energy}, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() } }); - case CALC_ES: - return TextFormat.sprintf("%3d %3d %-5s %-5s %6.3f %8.3f %8.3f", - new Object[] { minAtoms[c.iData[0]].sType, minAtoms[c.iData[1]].sType, - new float[] { (float)c.rab, (float)c.dData[0]/*qq*/, (float)c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() } }); + default: + return super.getDebugLine(iType, c); } - return ""; } - @Override - String getDebugFooter(int iType, double energy) { - String s = ""; - switch (iType){ - case CALC_DISTANCE: - s = "BOND STRETCHING"; - break; - case CALC_ANGLE: - s = "ANGLE BENDING"; - break; - case CALC_TORSION: - s = "TORSIONAL"; - break; - case CALC_OOP: - s = "OUT-OF-PLANE BENDING"; - break; - case CALC_VDW: - s = "VAN DER WAALS"; - break; - case CALC_ES: - s = "ELECTROSTATIC ENERGY"; - break; - } - return TextFormat.sprintf("\n TOTAL %s ENERGY = %8.3f %s\n", - new Object[] { s, getUnit(), Float.valueOf((float) energy) }); - } } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java 2012-05-13 04:29:32 UTC (rev 17125) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java 2012-05-14 03:17:17 UTC (rev 17126) @@ -33,7 +33,6 @@ import org.jmol.minimize.MinBond; import org.jmol.minimize.MinTorsion; import org.jmol.minimize.Util; -import org.jmol.util.TextFormat; /* * Java implementation by Bob Hanson 3/2008 @@ -142,7 +141,7 @@ return true; } - private boolean isInvertible(int n) { + private static boolean isInvertible(int n) { switch (n) { case 6: // C case 7: // N @@ -218,17 +217,8 @@ double compute(Object[] dataIn) { getPointers(dataIn); r0 = dData[0]; - kb = dData[1]; - ia = iData[0]; - ib = iData[1]; - - if (gradients) { - da.set(minAtoms[ia].coord); - db.set(minAtoms[ib].coord); - rab = Util.restorativeForceAndDistance(da, db, dc); - } else { - rab = Math.sqrt(Util.distance2(minAtoms[ia].coord, minAtoms[ib].coord)); - } + kb = dData[1]; + setBondVariables(this); // Er = 0.5 k (r - r0)^2 @@ -237,8 +227,7 @@ if (gradients) { dE = 2.0 * kb * delta; - addForce(da, ia, dE); - addForce(db, ib, dE); + addForces(this, 2); } if (logging) @@ -308,34 +297,20 @@ * (3.0 * rab * rbc * (1.0 - cosT0 * cosT0) - rac * rac * cosT0); calc.add(new Object[] { new int[] { ia, ib, ic, coordination }, - new double[] { ka, c0 - c2, c1, 2 * c2, theta0 * RAD_TO_DEG, preliminaryMagnification * ka } }); + new double[] { ka, theta0 * RAD_TO_DEG, c0 - c2, c1, 2 * c2, preliminaryMagnification * ka } }); } @Override double compute(Object[] dataIn) { getPointers(dataIn); - ia = iData[0]; - ib = iData[1]; - ic = iData[2]; int coordination = iData[3]; - double ka = (isPreliminary ? dData[4] : dData[0]); - double a0 = dData[1]; - double a1 = dData[2]; - double a2 = dData[3]; - - if (gradients) { - da.set(minAtoms[ia].coord); - db.set(minAtoms[ib].coord); - dc.set(minAtoms[ic].coord); - theta = Util.restorativeForceAndAngleRadians(da, db, dc); - } else { - theta = Util.getAngleRadiansABC(minAtoms[ia].coord, minAtoms[ib].coord, minAtoms[ic].coord); - } + double ka = (isPreliminary ? dData[5] : dData[0]); + double a0 = dData[2]; + double a1 = dData[3]; + double a2 = dData[4]; + setAngleVariables(this); - if (!Util.isFinite(theta)) - theta = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us; - //problem here for square planar cis or trans if ((coordination == 4 || coordination == 6) && (theta > 2.35619 || theta < 0.785398)) // 135o, 45o @@ -377,9 +352,7 @@ default: dE = -ka * (a1 * sinT - 2.0 * a2 * cosT * sinT); } - addForce(da, ia, dE); - addForce(db, ib, dE); - addForce(dc, ic, dE); + addForces(this, 3); } if (logging) @@ -496,36 +469,16 @@ double compute(Object[] dataIn) { getPointers(dataIn); - - double V = dData[0]; - double cosNPhi0 = dData[1]; - ia = iData[0]; - ib = iData[1]; - ic = iData[2]; - id = iData[3]; int n = iData[4]; - - if (gradients) { - da.set(minAtoms[ia].coord); - db.set(minAtoms[ib].coord); - dc.set(minAtoms[ic].coord); - dd.set(minAtoms[id].coord); - theta = Util.restorativeForceAndTorsionAngleRadians(da, db, dc, dd); - if (!Util.isFinite(theta)) - theta = 0.001 * DEG_TO_RAD; - } else { - theta = Util.getTorsionAngleRadians(minAtoms[ia].coord, minAtoms[ib].coord, - minAtoms[ic].coord, minAtoms[id].coord, v1, v2, v3); - } + double V = dData[0]; + double cosNPhi0 = dData[1]; + setTorsionVariables(this); energy = V * (1.0 - cosNPhi0 * Math.cos(theta * n)); if (gradients) { dE = V * n * cosNPhi0 * Math.sin(n * theta); - addForce(da, ia, dE); - addForce(db, ib, dE); - addForce(dc, ic, dE); - addForce(dd, id, dE); + addForces(this, 4); } if (logging) @@ -733,33 +686,14 @@ double a0 = dData[1]; double a1 = dData[2]; double a2 = dData[3]; - ia = iData[0]; - ib = iData[1]; - ic = iData[2]; - id = iData[3]; + setOopVariables(this); + + double cosTheta = Math.cos(theta); - da.set(minAtoms[ia].coord); - db.set(minAtoms[ib].coord); - dc.set(minAtoms[ic].coord); - dd.set(minAtoms[id].coord); - - if (gradients) { - theta = Util.restorativeForceAndOutOfPlaneAngleRadians(da, db, dc, dd, v1, v2, v3); - } else { - //if (atoms[ib].atom.getAtomIndex() == 20) - //System.out.println(" atom palne angl " + atoms[ib].atom.getInfo()); - - theta = Util.pointPlaneAngleRadians(da, db, dc, dd, v1, v2, v3); - } - - if (!Util.isFinite(theta)) - theta = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us; - //energy = koop * (c0 + c1 * Math.cos(theta) + c2 * Math.cos(2.0 * theta)); // //using - double cosTheta = Math.cos(theta); energy = koop * (a0 + a1 * cosTheta + a2 * cosTheta * cosTheta); //if (atoms[ib].atom.getAtomIndex() == 20) @@ -770,10 +704,7 @@ // not checked in Java dE = koop * (a1 * Math.sin(theta) + a2 * 2.0 * Math.sin(theta) * cosTheta); - addForce(da, ia, dE); - addForce(db, ib, dE); - addForce(dc, ic, dE); - addForce(dd, id, dE); + addForces(this, 4); } if (logging) @@ -821,19 +752,8 @@ getPointers(dataIn); double Xab = dData[0]; double Dab = dData[1]; - ia = iData[0]; - ib = iData[1]; - if (gradients) { - da.set(minAtoms[ia].coord); - db.set(minAtoms[ib].coord); - rab = Util.restorativeForceAndDistance(da, db, dc); - } else - rab = Math.sqrt(Util.distance2(minAtoms[ia].coord, minAtoms[ib].coord)); - - if (Util.isNearZero(rab, 1.0e-3)) - rab = 1.0e-3; - + setPairVariables(this); // Evdw = Dab [(Xab/r)^12 - 2(Xab/r)^6] Lennard-Jones // = Dab (Xab/r)^6[(Xab/r)^6 - 2] @@ -845,8 +765,7 @@ if (gradients) { dE = Dab * 12.0 * (1.0 - term6) * term6 * term / Xab; // unchecked - addForce(da, ia, dE); - addForce(db, ib, dE); + addForces(this, 2); } if (logging) @@ -873,29 +792,16 @@ } @Override - double compute(Object[] dataIn) { - + double compute(Object[] dataIn) { getPointers(dataIn); - double qq = dData[0]; - ia = iData[0]; - ib = iData[1]; - - if (gradients) { - da.set(minAtoms[ia].coord); - db.set(minAtoms[ib].coord); - rab = Util.restorativeForceAndDistance(da, db, dc); - } else - rab = Math.sqrt(Util.distance2(minAtoms[ia].coord, minAtoms[ib].coord)); + double qq = dData[0]; + setPairVariables(this); - if (Util.isNearZero(rab, 1.0e-3)) - rab = 1.0e-3; - energy = qq / rab; if (gradients) { dE = -qq / (rab * rab); - addForce(da, ia, dE); - addForce(db, ib, dE); + addForces(this, 2); } if (logging) @@ -909,158 +815,15 @@ ///////// REPORTING ///////////// @Override - String getAtomList(String title) { - String trailer = - "----------------------------------------" - + "-------------------------------------------------------\n"; - StringBuffer sb = new StringBuffer(); - sb.append("\n" + title + "\n\n" - + " ATOM X Y Z TYPE GRADX GRADY GRADZ " - + "---------BONDED ATOMS--------\n" - + trailer); - for (int i = 0; i < atomCount; i++) { - MinAtom atom = minAtoms[i]; - int[] others = atom.getBondedAtomIndexes(); - int[] iVal = new int[others.length + 1]; - iVal[0] = atom.atom.getAtomNumber(); - String s = " "; - for (int j = 0; j < others.length; j++) { - s += " %3d"; - iVal[j + 1] = minAtoms[others[j]].atom.getAtomNumber(); - } - sb.append(TextFormat.sprintf("%3d %8.3f %8.3f %8.3f %-5s %8.3f %8.3f %8.3f" + s + "\n", - new Object[] { atom.sType, - new float[] { (float) atom.coord[0], (float) atom.coord[1], - (float) atom.coord[2], (float) atom.force[0], (float) atom.force[1], - (float) atom.force[2] }, iVal})); - } - sb.append(trailer + "\n\n"); - return sb.toString(); - } - - @Override String getDebugHeader(int iType) { switch (iType){ case -1: return "Universal Force Field -- " + "Rappe, A. K., et. al.; J. Am. Chem. Soc. (1992) 114(25) p. 10024-10035\n"; - case CALC_DISTANCE: - return - "\nB O N D S T R E T C H I N G (" + bondCount + " bonds)\n\n" - +" ATOMS ATOM TYPES BOND BOND IDEAL FORCE\n" - +" I J I J TYPE LENGTH LENGTH CONSTANT DELTA ENERGY\n" - +"--------------------------------------------------------------------------------"; - case CALC_ANGLE: - return - "\nA N G L E B E N D I N G (" + minAngles.length + " angles)\n\n" - +" ATOMS ATOM TYPES VALENCE IDEAL FORCE\n" - +" I J K I J K ANGLE ANGLE CONSTANT ENERGY\n" - +"--------------------------------------------------------------------------"; - case CALC_TORSION: - return - "\nT O R S I O N A L (" + minTorsions.length + " torsions)\n\n" - +" ATOMS ATOM TYPES n COS FORCE TORSION\n" - +" I J K L I J K L (n phi0) CONSTANT ANGLE ENERGY\n" - +"---------------------------------------------------------------------------------------------"; - case CALC_OOP: - return - "\nO U T - O F - P L A N E B E N D I N G\n\n" - +" ATOMS ATOM TYPES OOP FORCE \n" - +" I J K L I J K L ANGLE CONSTANT ENERGY\n" - +"--------------------------------------------------------------------------"; - case CALC_VDW: - return - "\nV A N D E R W A A L S\n\n" - +" ATOMS ATOM TYPES\n" - +" I J I J Rij kij ENERGY\n" - +"-----------------------------------------------"; - case CALC_ES: - return - "\nE L E C T R O S T A T I C I N T E R A C T I O N S\n\n" - +" ATOMS ATOM TYPES QiQj\n" - +" I J I J Rij *332.17 ENERGY\n" - +"-----------------------------------------------"; + default: + return super.getDebugHeader(iType); } - return ""; } - String getDebugLine(int iType, Calculation c) { - switch (iType) { - case CALC_DISTANCE: - return TextFormat.sprintf( - "%3d %3d %-5s %-5s %4.2f%8.3f %8.3f %8.3f %8.3f %8.3f", - new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, - new float[] { (float)c.dData[2]/*rab*/, (float)c.rab, - (float)c.dData[0], (float)c.dData[1], - (float)c.delta, (float)c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() }}); - case CALC_ANGLE: - return TextFormat.sprintf( - "%3d %3d %3d %-5s %-5s %-5s %8.3f %8.3f %8.3f %8.3f", - new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, - minAtoms[c.ic].sType, - new float[] { (float)(c.theta * RAD_TO_DEG), (float)c.dData[4] /*THETA0*/, - (float)c.dData[0]/*Kijk*/, (float) c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), - minAtoms[c.ic].atom.getAtomNumber()} }); - case CALC_TORSION: - return TextFormat.sprintf( - "%3d %3d %3d %3d %-5s %-5s %-5s %-5s %3d %8.3f %8.3f %8.3f %8.3f", - new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, - minAtoms[c.ic].sType, minAtoms[c.id].sType, - new float[] { (float) c.dData[1]/*cosNphi0*/, (float) c.dData[0]/*V*/, - (float) (c.theta * RAD_TO_DEG), (float) c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), - minAtoms[c.ic].atom.getAtomNumber(), minAtoms[c.id].atom.getAtomNumber(), c.iData[4] } }); - case CALC_OOP: - return TextFormat.sprintf("" + - "%3d %3d %3d %3d %-5s %-5s %-5s %-5s %8.3f %8.3f %8.3f", - new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, - minAtoms[c.ic].sType, minAtoms[c.id].sType, - new float[] { (float)(c.theta * RAD_TO_DEG), - (float)c.dData[0]/*koop*/, (float) c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber(), - minAtoms[c.ic].atom.getAtomNumber(), minAtoms[c.id].atom.getAtomNumber() } }); - case CALC_VDW: - return TextFormat.sprintf("%3d %3d %-5s %-5s %6.3f %8.3f %8.3f", - new Object[] { minAtoms[c.iData[0]].sType, minAtoms[c.iData[1]].sType, - new float[] { (float)c.rab, (float)c.dData[0]/*kab*/, (float)c.energy}, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() } }); - case CALC_ES: - return TextFormat.sprintf("%3d %3d %-5s %-5s %6.3f %8.3f %8.3f", - new Object[] { minAtoms[c.iData[0]].sType, minAtoms[c.iData[1]].sType, - new float[] { (float)c.rab, (float)c.dData[0]/*qq*/, (float)c.energy }, - new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() } }); - } - return ""; - } - - @Override - String getDebugFooter(int iType, double energy) { - String s = ""; - switch (iType){ - case CALC_DISTANCE: - s = "BOND STRETCHING"; - break; - case CALC_ANGLE: - s = "ANGLE BENDING"; - break; - case CALC_TORSION: - s = "TORSIONAL"; - break; - case CALC_OOP: - s = "OUT-OF-PLANE BENDING"; - break; - case CALC_VDW: - s = "VAN DER WAALS"; - break; - case CALC_ES: - s = "ELECTROSTATIC ENERGY"; - break; - } - return TextFormat.sprintf("\n TOTAL %s ENERGY = %8.3f %s\n", - new Object[] { s, getUnit(), Float.valueOf((float) energy) }); - } - } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java 2012-05-13 04:29:32 UTC (rev 17125) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java 2012-05-14 03:17:17 UTC (rev 17126) @@ -578,7 +578,8 @@ for (int i = minAtomCount; --i >= 0;) { int it = rawAtomTypes[minAtoms[i].atom.index]; minAtoms[i].ffAtomType = atomTypes.get(Math.max(0, it)); - minAtoms[i].ffType = (it < 0 ? -it : atomTypes.get(it).mmType); + minAtoms[i].ffType = (it < 0 ? -it : atomTypes.get(it).mmType); + minAtoms[i].vdwKey = MinObject.getKey(KEY_VDW, minAtoms[i].ffType, 127, 127, A4_VDW); } // fix order in bonds and set type and key @@ -783,45 +784,11 @@ } double getOutOfPlaneParameter(int[] data) { - // defined for typeB = 2, 3, 8, 10, 17, 26, 30, 37, 39, 40, 41, 43, - // 45, 49, 54, 55, 56, 57, 58, 63, 64, 67, 69, 73, 78, 80, 81, 82 - // but is 0 for 8 (amines), 17 (sulfoxides), - // 26 (PD3), 43 (N-S), 73 (O-S(=O)R, 82 (N-oxide) - // that is, just the planar ones: - // 2, 3, 10, 30, 37, 39, 40, 41, - // 45, 49, 54, 55, 56, 57, 58, 63, - // 64, 67, 69, 78, 80, 81 - int typeB = minAtoms[data[1]].ffType; - switch (typeB) { - default: - return 0; - case 2: - case 3: - case 10: - case 30: - case 37: - case 39: - case 40: - case 41: - case 45: - case 49: - case 54: - case 55: - case 56: - case 57: - case 58: - case 63: - case 64: - case 67: - case 69: - case 78: - case 80: - case 81: - break; - } fixOrder(data, 0, 2); + fixOrder(data, 0, 3); fixOrder(data, 2, 3); int typeA = minAtoms[data[0]].ffType; + int typeB = minAtoms[data[1]].ffType; int typeC = minAtoms[data[2]].ffType; int typeD = minAtoms[data[3]].ffType; Object params = getFFParams(MinObject.getKey(KEY_OOP, This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2012-05-14 22:39:02
|
Revision: 17129 http://jmol.svn.sourceforge.net/jmol/?rev=17129&view=rev Author: hansonr Date: 2012-05-14 22:38:55 +0000 (Mon, 14 May 2012) Log Message: ----------- MMFF -- preliminary parameter setting only for -- all except electrostatic now set forceField "MMFF" minimize energy Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/MinAngle.java trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java Modified: trunk/Jmol/src/org/jmol/minimize/MinAngle.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/MinAngle.java 2012-05-14 15:19:52 UTC (rev 17128) +++ trunk/Jmol/src/org/jmol/minimize/MinAngle.java 2012-05-14 22:38:55 UTC (rev 17129) @@ -4,6 +4,6 @@ public int sbType; public Integer sbKey; MinAngle(int[] data) { - this.data = data; + this.data = data; // ia, ib, ic, iab, ibc } } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java 2012-05-14 15:19:52 UTC (rev 17128) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java 2012-05-14 22:38:55 UTC (rev 17129) @@ -125,6 +125,7 @@ } void appendLogData(String s) { + System.out.println("calcultions LOG: " + s); logData.append(s).append("\n"); } @@ -447,6 +448,12 @@ +" ATOMS ATOM TYPES VALENCE IDEAL FORCE\n" +" I J K I J K ANGLE ANGLE CONSTANT ENERGY\n" +"--------------------------------------------------------------------------"; + case CALC_STRETCH_BEND: + return + "\nS T R E T C H B E N D I N G (" + (minAngles.length * 2) + " angles)\n\n" + +" ATOMS ATOM TYPES VALENCE IDEAL FORCE\n" + +" I J K I J K ANGLE ANGLE CONSTANT ENERGY\n" + +"--------------------------------------------------------------------------"; case CALC_TORSION: return "\nT O R S I O N A L (" + minTorsions.length + " torsions)\n\n" @@ -486,6 +493,7 @@ (float)c.delta, (float)c.energy }, new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() }}); case CALC_ANGLE: + case CALC_STRETCH_BEND: return TextFormat.sprintf( "%3d %3d %3d %-5s %-5s %-5s %8.3f %8.3f %8.3f %8.3f", new Object[] { minAtoms[c.ia].sType, minAtoms[c.ib].sType, @@ -611,7 +619,7 @@ c.rab = 1.0e-3; } - private void setCoords(Calculation c, int n) { + void setCoords(Calculation c, int n) { switch(n) { case 4: da.set(minAtoms[c.ia].coord); Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java 2012-05-14 15:19:52 UTC (rev 17128) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java 2012-05-14 22:38:55 UTC (rev 17129) @@ -294,12 +294,53 @@ class SBCalc extends Calculation { - // TODO void setData(List<Object[]> calc, MinAngle angle) { + // not applicable for linear types + switch (minAtoms[angle.data[1]].ffType) { + case 4: + case 53: + case 61: + return; + } + double[] data = (double[]) getParameter(angle.sbKey); + double[] datakat0 = (double[]) getParameter(angle.key); + double[] dataij = (double[]) getParameter(minBonds[angle.data[3]].key); + double[] datajk = (double[]) getParameter(minBonds[angle.data[4]].key); + if (data == null || datakat0 == null || dataij == null || datajk == null) + return; + double theta0 = datakat0[1]; + double r0ij = dataij[1]; + double r0jk = datajk[1]; + calc.add(new Object[] { angle.data, new double[] { data[0], theta0, r0ij } }); + calc.add(new Object[] { new int[] {angle.data[2], angle.data[1], angle.data[0]}, + new double[] { data[1], theta0, r0jk } }); } @Override double compute(Object[] dataIn) { + getPointers(dataIn); + double k = 2.51210 * dData[0]; + double t0 = dData[1]; + double r0_ab = dData[2]; + + setBondVariables(this); + setAngleVariables(this); + double dr_ab = rab - r0_ab; + delta = theta * RAD_TO_DEG - t0; + // equation 5 + energy = k * dr_ab * delta; + + if (logging) + appendLogData(getDebugLine(CALC_STRETCH_BEND, this)); + + if (gradients) { + dE = k * dr_ab; + addForces(this, 3); + setBondVariables(this); + dE = k * delta; + addForces(this, 2); + } + return energy; } } @@ -412,18 +453,19 @@ double[] dataB = (double[]) getParameter(b.vdwKey); if (dataA == null || dataB == null) return; + double alpha_a = dataA[0]; double N_a = dataA[1]; double A_a = dataA[2]; double G_a = dataA[3]; - int DA_a = (int) dataA[4]; + int DA_a = (int) dataA[4]; + double alpha_b = dataB[0]; double N_b = dataB[1]; double A_b = dataB[2]; double G_b = dataB[3]; int DA_b = (int) dataB[4]; - double rs_aa = A_a * Math.pow(alpha_a, 0.25); double rs_bb = A_b * Math.pow(alpha_b, 0.25); double gamma = (rs_aa - rs_bb) / (rs_aa + rs_bb); Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java 2012-05-14 15:19:52 UTC (rev 17128) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java 2012-05-14 22:38:55 UTC (rev 17129) @@ -50,8 +50,8 @@ public class ForceFieldMMFF extends ForceField { - private static final int A4_SB = 1; - private static final int A4_SBDEF = 2; + private static final int A4_SB = 125; + private static final int A4_SBDEF = 126; private static final int A4_VDW = 3; private static final int A4_CHRG = 4; private static final int KEY_VDW = 0; @@ -72,7 +72,6 @@ return rawMMFF94Charges; } - /* * from SMARTS search when calculating partial charges: * @@ -681,7 +680,7 @@ private int getRowFor(int i) { int elemno = minAtoms[i].atom.getElementNumber(); - return (elemno < 3 ? 1 : elemno < 11 ? 2 : elemno < 19 ? 3 : 4); + return (elemno < 3 ? 0 : elemno < 11 ? 1 : elemno < 19 ? 2 : elemno < 37 ? 3 : 4); } private static Object getFFParams(Integer key) { This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2012-05-15 04:27:47
|
Revision: 17134 http://jmol.svn.sourceforge.net/jmol/?rev=17134&view=rev Author: hansonr Date: 2012-05-15 04:27:40 +0000 (Tue, 15 May 2012) Log Message: ----------- MMFF94 full single-point energy calculation. Minimization itself not tested; reasonably good validation; result is extremely sensitive to electrostatic charge; not sure what is exactly going on there. Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/MinAtom.java trunk/Jmol/src/org/jmol/minimize/Minimizer.java trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldUFF.java Modified: trunk/Jmol/src/org/jmol/minimize/MinAtom.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/MinAtom.java 2012-05-14 23:56:19 UTC (rev 17133) +++ trunk/Jmol/src/org/jmol/minimize/MinAtom.java 2012-05-15 04:27:40 UTC (rev 17134) @@ -25,6 +25,7 @@ package org.jmol.minimize; import java.util.ArrayList; +import java.util.BitSet; import java.util.List; import org.jmol.minimize.forcefield.AtomType; @@ -32,7 +33,7 @@ public class MinAtom { - int rawIndex; + int index; public String sType; public Atom atom; public AtomType ffAtomType; @@ -40,22 +41,26 @@ public Integer vdwKey; public double[] coord = new double[3]; public double[] force = new double[3]; - public List<MinBond> bonds = new ArrayList<MinBond>(); + private List<MinBond> bonds = new ArrayList<MinBond>(); public int nBonds; public int hCount; public double partialCharge; + public BitSet bsVdw = new BitSet(); + public BitSet bs14 = new BitSet(); int[] bondedAtoms; @Override public String toString() { - return "#" + rawIndex + " " + sType; + return "#" + index + " " + sType; } - MinAtom(int index, Atom atom, double[] coord) { - this.rawIndex = index; + MinAtom(int index, Atom atom, double[] coord, int atomCount) { + this.index = index; this.atom = atom; this.coord = coord; + bsVdw.set(index + 1, atomCount); + bsVdw.clear(index); hCount = atom.getCovalentHydrogenCount(); } @@ -77,7 +82,7 @@ if (bondedAtoms == null) { bondedAtoms = new int[nBonds]; for (int i = nBonds; --i >= 0;) - bondedAtoms[i] = bonds.get(i).getOtherAtom(rawIndex); + bondedAtoms[i] = bonds.get(i).getOtherAtom(index); } return bondedAtoms; } @@ -86,4 +91,14 @@ return atom.getInfo(); } + public void addBond(MinBond bond, int i) { + bonds.add(bond); + nBonds++; + bsVdw.clear(i); + } + + public int getBondIndex(int j) { + return bonds.get(j).index; + } + } Modified: trunk/Jmol/src/org/jmol/minimize/Minimizer.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2012-05-14 23:56:19 UTC (rev 17133) +++ trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2012-05-15 04:27:40 UTC (rev 17134) @@ -287,7 +287,7 @@ elemnoMax = Math.max(elemnoMax, atomicNo); bsElements.set(atomicNo); minAtoms[pt] = new MinAtom(pt, atom, new double[] { atom.x, atom.y, - atom.z }); + atom.z }, atomCount); minAtoms[pt].sType = atom.getAtomName(); } @@ -353,10 +353,8 @@ MinBond bond = minBonds[i] = bondInfo.get(i); int atom1 = bond.data[0]; int atom2 = bond.data[1]; - minAtoms[atom1].bonds.add(bond); - minAtoms[atom2].bonds.add(bond); - minAtoms[atom1].nBonds++; - minAtoms[atom2].nBonds++; + minAtoms[atom1].addBond(bond, atom2); + minAtoms[atom2].addBond(bond, atom1); } for (int i = 0; i < atomCount; i++) minAtoms[i].getBondedAtomIndexes(); @@ -375,7 +373,8 @@ for (int j = atomList.length; --j >= 0;) if ((ic = atomList[j]) > ia) { vAngles.add(new MinAngle(new int[] { ia, ib, ic, i, - minAtoms[ib].bonds.get(j).index})); + minAtoms[ib].getBondIndex(j)})); + minAtoms[ia].bsVdw.clear(ic); /* System.out.println (" " + minAtoms[ia].getIdentity() + " -- " + minAtoms[ib].getIdentity() + " -- " @@ -388,8 +387,9 @@ for (int j = atomList.length; --j >= 0;) if ((ic = atomList[j]) < ib && ic > ia) { vAngles - .add(new MinAngle(new int[] { ic, ia, ib, minAtoms[ia].bonds.get(j).index, + .add(new MinAngle(new int[] { ic, ia, ib, minAtoms[ia].getBondIndex(j), i})); + minAtoms[ic].bsVdw.clear(ib); System.out.println ("a " + minAtoms[ic].getIdentity() + " -- " + minAtoms[ia].getIdentity() + " -- " @@ -419,7 +419,8 @@ if (id != ia && id != ib) { vTorsions.add(new MinTorsion(new int[] { ia, ib, ic, id, angle[ForceField.ABI_IJ], angle[ForceField.ABI_JK], - minAtoms[ic].bonds.get(j).index })); + minAtoms[ic].getBondIndex(j) })); + minAtoms[Math.min(ia, id)].bs14.set(Math.max(ia, id)); /* System.out.println("t " + minAtoms[ia].getIdentity() + " -- " + minAtoms[ib].getIdentity() + " -- " + minAtoms[ic].getIdentity() + " -- " @@ -435,7 +436,8 @@ if (id != ic && id != ib) { vTorsions.add(new MinTorsion(new int[] { ic, ib, ia, id, angle[ForceField.ABI_JK], angle[ForceField.ABI_IJ], - minAtoms[ia].bonds.get(j).index })); + minAtoms[ia].getBondIndex(j) })); + minAtoms[Math.min(ic, id)].bs14.set(Math.max(ic, id)); /* System.out.println("t " + minAtoms[ic].getIdentity() + " -- " + minAtoms[ib].getIdentity() + " -- " + minAtoms[ia].getIdentity() + " -- " Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java 2012-05-14 23:56:19 UTC (rev 17133) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/Calculations.java 2012-05-15 04:27:40 UTC (rev 17134) @@ -24,6 +24,7 @@ package org.jmol.minimize.forcefield; +import java.util.BitSet; import java.util.List; import javax.vecmath.Vector3d; @@ -63,8 +64,6 @@ MinBond[] minBonds; MinAngle[] minAngles; MinTorsion[] minTorsions; - double[] partialCharges; - boolean havePartialCharges; List<Object[]> constraints; boolean isPreliminary; @@ -75,7 +74,6 @@ Calculations(ForceField ff, MinAtom[] minAtoms, MinBond[] minBonds, MinAngle[] minAngles, MinTorsion[] minTorsions, - double[] partialCharges, List<Object[]> constraints) { this.ff = ff; this.minAtoms = minAtoms; @@ -87,15 +85,6 @@ angleCount = minAngles.length; torsionCount = minTorsions.length; this.constraints = constraints; - if (partialCharges != null && partialCharges.length == atomCount) - for (int i = atomCount; --i >= 0;) - if (partialCharges[i] != 0) { - havePartialCharges = true; - break; - } - if (!havePartialCharges) - partialCharges = null; - this.partialCharges = partialCharges; } abstract boolean setupCalculations(); @@ -125,7 +114,6 @@ } void appendLogData(String s) { - System.out.println("calcultions LOG: " + s); logData.append(s).append("\n"); } @@ -150,37 +138,9 @@ protected void pairSearch(List<Object[]> calc1, PairCalc pc1, List<Object[]> calc2, PairCalc pc2) { - /*A:*/ for (int i = 0; i < atomCount - 1; i++) { // one atom... - MinAtom atomA = minAtoms[i]; - int[] atomList1 = atomA.getBondedAtomIndexes(); - B: for (int j = i + 1; j < atomCount; j++) { // another atom... - MinAtom atomB = minAtoms[j]; - /*nbrA:*/ for (int k = atomList1.length; --k >= 0;) { // check bonded A-B - MinAtom nbrA = minAtoms[atomList1[k]]; - if (nbrA == atomB) - continue B; // pick another B - if (nbrA.nBonds == 1) - continue; - int[] atomList2 = nbrA.getBondedAtomIndexes(); // check A-X-B - /*nbrAA:*/ for (int l = atomList2.length; --l >= 0;) { - MinAtom nbrAA = minAtoms[atomList2[l]]; - if (nbrAA == atomB) - continue B; // pick another B - - //this next would exclude A-X-X-B, but Rappe does not do that, he says - -/* if (nbrAA.nBonds == 1) - continue; - int[] atomList3 = nbrAA.getBondedAtomIndexes(); // check A-X-X-B - nbrAAA: for (int m = atomList3.length; --m >= 0;) { - MinAtom nbrAAA = atoms[atomList3[m]]; - if (nbrAAA == atomB) - continue B; // pick another B - } -*/ - - } - } + for (int i = 0; i < atomCount - 1; i++) { + BitSet bsVdw = minAtoms[i].bsVdw; + for (int j = bsVdw.nextSetBit(0); j >= 0; j = bsVdw.nextSetBit(j + 1)) { pc1.setData(calc1, i, j); if (pc2 != null) pc2.setData(calc2, i, j); @@ -475,9 +435,9 @@ case CALC_ES: return "\nE L E C T R O S T A T I C I N T E R A C T I O N S\n\n" - +" ATOMS ATOM TYPES QiQj\n" - +" I J I J Rij *332.17 ENERGY\n" - +"-----------------------------------------------"; + +" ATOMS ATOM TYPES \n" + +" I J I J Rij f Qi Qj ENERGY\n" + +"-------------------------------------------------------------------"; } return ""; } @@ -526,9 +486,9 @@ new float[] { (float)c.rab, (float)c.dData[0]/*kab*/, (float)c.energy}, new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() } }); case CALC_ES: - return TextFormat.sprintf("%3d %3d %-5s %-5s %6.3f %8.3f %8.3f", + return TextFormat.sprintf("%3d %3d %-5s %-5s %6.3f %8.3f %8.3f %8.3f %8.3f", new Object[] { minAtoms[c.iData[0]].sType, minAtoms[c.iData[1]].sType, - new float[] { (float)c.rab, (float)c.dData[0]/*qq*/, (float)c.energy }, + new float[] { (float)c.rab, (float)c.dData[0]/*q1*/, (float)c.dData[1]/*q2*/, (float)c.dData[2]/*f*/, (float)c.energy }, new int[] { minAtoms[c.ia].atom.getAtomNumber(), minAtoms[c.ib].atom.getAtomNumber() } }); } return ""; @@ -562,15 +522,17 @@ //////////////////////////////////////////////////// - void setBondVariables(Calculation c) { + void setPairVariables(Calculation c) { if (gradients) { setCoords(c, 2); c.rab = Util.restorativeForceAndDistance(da, db, dc); } else { c.rab = Math.sqrt(Util.distance2(minAtoms[c.ia].coord, minAtoms[c.ib].coord)); } + if (Util.isNearZero(c.rab, 1.0e-3)) + c.rab = 1.0e-3; } - + void setAngleVariables(Calculation c) { if (gradients) { setCoords(c, 3); @@ -587,9 +549,6 @@ if (gradients) { c.theta = Util.restorativeForceAndOutOfPlaneAngleRadians(da, db, dc, dd, v1, v2, v3); } else { - //if (atoms[ib].atom.getAtomIndex() == 20) - //System.out.println(" atom palne angl " + atoms[ib].atom.getInfo()); - c.theta = Util.pointPlaneAngleRadians(da, db, dc, dd, v1, v2, v3); } if (!Util.isFinite(c.theta)) @@ -608,17 +567,6 @@ } } - void setPairVariables(Calculation c) { - if (gradients) { - setCoords(c, 2); - c.rab = Util.restorativeForceAndDistance(da, db, dc); - } else { - c.rab = Math.sqrt(Util.distance2(minAtoms[c.ia].coord, minAtoms[c.ib].coord)); - } - if (Util.isNearZero(c.rab, 1.0e-3)) - c.rab = 1.0e-3; - } - void setCoords(Calculation c, int n) { switch(n) { case 4: Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java 2012-05-14 23:56:19 UTC (rev 17133) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java 2012-05-15 04:27:40 UTC (rev 17134) @@ -32,7 +32,6 @@ import org.jmol.minimize.MinAtom; import org.jmol.minimize.MinBond; import org.jmol.minimize.MinTorsion; -import org.jmol.minimize.Util; import org.jmol.util.TextFormat; /* @@ -82,9 +81,9 @@ CalculationsMMFF(ForceField ff, Map<Integer, Object> ffParams, MinAtom[] minAtoms, MinBond[] minBonds, - MinAngle[] minAngles, MinTorsion[] minTorsions, double[] partialCharges, + MinAngle[] minAngles, MinTorsion[] minTorsions, List<Object[]> constraints) { - super(ff, minAtoms, minBonds, minAngles, minTorsions, partialCharges, constraints); + super(ff, minAtoms, minBonds, minAngles, minTorsions, constraints); this.ffParams = ffParams; bondCalc = new DistanceCalc(); angleCalc = new AngleCalc(); @@ -228,7 +227,7 @@ getPointers(dataIn); kb = dData[0]; r0 = dData[1]; - setBondVariables(this); + setPairVariables(this); delta = rab - r0; delta2 = delta * delta; @@ -323,7 +322,7 @@ double t0 = dData[1]; double r0_ab = dData[2]; - setBondVariables(this); + setPairVariables(this); setAngleVariables(this); double dr_ab = rab - r0_ab; delta = theta * RAD_TO_DEG - t0; @@ -336,7 +335,7 @@ if (gradients) { dE = k * dr_ab; addForces(this, 3); - setBondVariables(this); + setPairVariables(this); dE = k * delta; addForces(this, 2); } @@ -485,10 +484,9 @@ @Override double compute(Object[] dataIn) { getPointers(dataIn); + setPairVariables(this); double rs = dData[0]; double eps = dData[1]; - - rab = Math.sqrt(Util.distance2(minAtoms[ia].coord, minAtoms[ib].coord)); double r_rs = rab / rs; energy = eps * Math.pow(1.07/ (r_rs + 0.07), 7) * (1.12 / (Math.pow(r_rs, 7) + 0.12) - 2); @@ -511,14 +509,34 @@ class ESCalc extends PairCalc { + private static final double BUFF = 0.05; + @Override void setData(List<Object[]> calc, int ia, int ib) { - // TODO -- use torsions for 1,4 business + if (minAtoms[ia].partialCharge == 0 || minAtoms[ib].partialCharge == 0) + return; + calc.add(new Object[] { new int[] { ia, ib }, new double[] { + minAtoms[ia].partialCharge, minAtoms[ib].partialCharge, + (minAtoms[ia].bs14.get(ib) ? 249.054 : 332.0716) } + }); } @Override double compute(Object[] dataIn) { - // TODO + getPointers(dataIn); + double f = dData[0] * dData[1] * dData[2]; + setPairVariables(this); + double d = rab + BUFF; + energy = f / d; // DIEL = 1 here + + if (gradients) { + dE = -energy / d; + addForces(this, 2); + } + + if (logging) + appendLogData(getDebugLine(CALC_ES, this)); + return energy; } } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java 2012-05-14 23:56:19 UTC (rev 17133) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsUFF.java 2012-05-15 04:27:40 UTC (rev 17134) @@ -71,20 +71,20 @@ TorsionCalc torsionCalc; OOPCalc oopCalc; VDWCalc vdwCalc; - ESCalc esCalc; + //ESCalc esCalc; CalculationsUFF(ForceField ff, Map<Object, FFParam> ffParams, MinAtom[] minAtoms, MinBond[] minBonds, MinAngle[] minAngles, MinTorsion[] minTorsions, - double[] partialCharges, List<Object[]> constraints) { - super(ff, minAtoms, minBonds, minAngles, minTorsions, partialCharges, constraints); + List<Object[]> constraints) { + super(ff, minAtoms, minBonds, minAngles, minTorsions, constraints); this.ffParams = ffParams; bondCalc = new DistanceCalc(); angleCalc = new AngleCalc(); torsionCalc = new TorsionCalc(); oopCalc = new OOPCalc(); vdwCalc = new VDWCalc(); - esCalc = new ESCalc(); + //esCalc = new ESCalc(); } @Override @@ -133,11 +133,7 @@ // it does not actually use it. Both Towhee and the UFF FAQ // discourage the use of electrostatics with UFF. - if (partialCharges == null) - pairSearch(calculations[CALC_VDW] = new ArrayList<Object[]>(), new VDWCalc(), null, null); - else - pairSearch(calculations[CALC_VDW] = new ArrayList<Object[]>(), new VDWCalc(), - calculations[CALC_ES] = new ArrayList<Object[]>(), new ESCalc()); + pairSearch(calculations[CALC_VDW] = new ArrayList<Object[]>(), new VDWCalc(), null, null); return true; } @@ -185,8 +181,8 @@ return oopCalc.compute(dataIn); case CALC_VDW: return vdwCalc.compute(dataIn); - case CALC_ES: - return esCalc.compute(dataIn); + //case CALC_ES: + //return esCalc.compute(dataIn); } return 0.0; } @@ -195,6 +191,8 @@ return ffParams.get(a); } + final static double KCAL332 = KCAL_TO_KJ * 332.0637; + class DistanceCalc extends Calculation { double r0, kb; @@ -218,7 +216,7 @@ getPointers(dataIn); r0 = dData[0]; kb = dData[1]; - setBondVariables(this); + setPairVariables(this); // Er = 0.5 k (r - r0)^2 @@ -696,9 +694,6 @@ energy = koop * (a0 + a1 * cosTheta + a2 * cosTheta * cosTheta); - //if (atoms[ib].atom.getAtomIndex() == 20) - //System.out.println(ib + " testing oop theta cosTheta=" + (theta * 180/Math.PI) + " " + cosTheta + " energy=" + energy + " koop=" + koop + " a0 a1 a2="+ a0 + " " + a1 + " "+ a2); - if (gradients) { // somehow we already get the -1 from the OOPDeriv -- so we'll omit it here // not checked in Java @@ -774,8 +769,7 @@ return energy; } } - - final static double KCAL332 = KCAL_TO_KJ * 332.0637; +/* class ESCalc extends PairCalc { @@ -811,7 +805,7 @@ } } - +*/ ///////// REPORTING ///////////// @Override Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java 2012-05-14 23:56:19 UTC (rev 17133) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java 2012-05-15 04:27:40 UTC (rev 17134) @@ -111,11 +111,8 @@ setArrays(m.atoms, m.bsAtoms, m.bonds, m.rawBondCount); setModelFields(); fixTypes(); - double[] partialCharges = new double[minAtoms.length]; - for (int i = m.bsAtoms.nextSetBit(0), pt = 0; i >= 0; i = m.bsAtoms.nextSetBit(i + 1)) - partialCharges[pt++] = rawMMFF94Charges[i]; calc = new CalculationsMMFF(this, ffParams, minAtoms, minBonds, - minAngles, minTorsions, partialCharges, minimizer.constraints); + minAngles, minTorsions, minimizer.constraints); calc.setLoggingEnabled(true); return calc.setupCalculations(); } @@ -551,8 +548,6 @@ Bond[] bonds = atoms[i].getBonds(); if (bonds != null) { types[i] = -atomTypes.get(types[bonds[0].getOtherAtom(atoms[i]).index]).hType; - //System.out.println(atoms[i] + " " + bonds[0].getOtherAtom(atoms[i]) + " " +atomTypes.get(types[bonds[0].getOtherAtom(atoms[i]).index]).mmType - // + " " + atomTypes.get(types[bonds[0].getOtherAtom(atoms[i]).index]).hType); } } if (Logger.debugging) @@ -585,10 +580,12 @@ private void fixTypes() { // set atom types in minAtoms for (int i = minAtomCount; --i >= 0;) { - int it = rawAtomTypes[minAtoms[i].atom.index]; + int rawIndex = minAtoms[i].atom.index; + int it = rawAtomTypes[rawIndex]; minAtoms[i].ffAtomType = atomTypes.get(Math.max(0, it)); minAtoms[i].ffType = (it < 0 ? -it : atomTypes.get(it).mmType); minAtoms[i].vdwKey = MinObject.getKey(KEY_VDW, minAtoms[i].ffType, 127, 127, A4_VDW); + minAtoms[i].partialCharge = rawMMFF94Charges[rawIndex]; } // fix order in bonds and set type and key Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldUFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldUFF.java 2012-05-14 23:56:19 UTC (rev 17133) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldUFF.java 2012-05-15 04:27:40 UTC (rev 17134) @@ -66,7 +66,7 @@ return false; setAtomTypes(bsElements, elemnoMax); calc = new CalculationsUFF(this, ffParams, minAtoms, minBonds, - minAngles, minTorsions, minimizer.partialCharges, minimizer.constraints); + minAngles, minTorsions, minimizer.constraints); return calc.setupCalculations(); } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2012-05-17 16:10:41
|
Revision: 17152 http://jmol.svn.sourceforge.net/jmol/?rev=17152&view=rev Author: hansonr Date: 2012-05-17 16:10:31 +0000 (Thu, 17 May 2012) Log Message: ----------- checking calculated energies for 761 models 1 COMKAQ E= -7.3250003 Eref= -7.6177 diff= 0.2926998 2 DUVHUX10 E= 64.759995 Eref= 64.082855 diff= 0.6771393 3 FORJIF E= 35.978 Eref= 35.833878 diff= 0.14412308 4 JADLIJ E= 25.104 Eref= 24.7038 diff= 0.4001999 5 KEPKIZ E= 61.192997 Eref= 61.816277 diff= 0.6232796 6 PHOSLA10 E= 111.232994 Eref= 112.07078 diff= 0.8377838 7 PHOSLB10 E= -93.479004 Eref= -92.64081 diff= 0.8381958 8 ERULE_01 E= -22.894 Eref= -21.515108 diff= 1.378891 9 ERULE_02 E= 29.357998 Eref= 29.799572 diff= 0.4415741 10 ERULE_03 E= -3.326 Eref= -2.9351802 diff= 0.3908198 11 ERULE_04 E= -2.596 Eref= -2.31007 diff= 0.28592992 12 ERULE_08 E= 33.444004 Eref= 34.41382 diff= 0.9698143 for 761 atoms, 12 have energies differences outside the range -0.1 to 0.1 with a standard deviation of 0.08684954 $ Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/Minimizer.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java trunk/Jmol/src/org/jmol/minimize/forcefield/mmff/validate/checkmm.spt Modified: trunk/Jmol/src/org/jmol/minimize/Minimizer.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2012-05-17 15:24:00 UTC (rev 17151) +++ trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2012-05-17 16:10:31 UTC (rev 17152) @@ -390,10 +390,12 @@ .add(new MinAngle(new int[] { ic, ia, ib, minAtoms[ia].getBondIndex(j), i})); minAtoms[ic].bsVdw.clear(ib); +/* System.out.println ("a " + minAtoms[ic].getIdentity() + " -- " + minAtoms[ia].getIdentity() + " -- " + minAtoms[ib].getIdentity()); +*/ } } } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java 2012-05-17 15:24:00 UTC (rev 17151) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceFieldMMFF.java 2012-05-17 16:10:31 UTC (rev 17152) @@ -146,12 +146,14 @@ private static final int TYPE_PBCI = 0x1; private static final int TYPE_VDW = 0x11; private static final int TYPE_CHRG = 0x22; - private static final int TYPE_BOND = 0x3; // 0011 - private static final int TYPE_ANGLE = 0x5; // 0101 - private static final int TYPE_SB = 0x15; // 1 0101 - private static final int TYPE_SBDEF = 0x25; // 10 0101 - private static final int TYPE_TORSION = 0x9; // 0 1001 - private static final int TYPE_OOP = 0xB; // 1 1011; + // the following are bit flags indicating in the 0xF range + // which atoms might have default parameter values + private static final int TYPE_BOND = 0x3; // 0011 + private static final int TYPE_ANGLE = 0x5; // 0101 + private static final int TYPE_SB = 0x15; // 01 0101 + private static final int TYPE_SBDEF = 0x25; // 10 0101 + private static final int TYPE_TORSION = 0x9; // 00 1001 + private static final int TYPE_OOP = 0xD; // 00 1101; private static List<AtomType> atomTypes; @@ -789,8 +791,9 @@ return false; } - private Integer getKey(MinObject o, int type, int ktype) { - int[] data = o.data; + private Integer getKey(Object obj, int type, int ktype) { + MinObject o = (obj instanceof MinObject ? (MinObject) obj : null); + int[] data = (o == null ? (int[]) obj : o.data); int n = 4; switch (ktype) { case TYPE_BOND: @@ -824,8 +827,14 @@ Integer key = null; for (int i = 0; i < 4; i++) typeData[i] = (i < n ? typeOf(data[i]) : 127); - if (ktype == TYPE_SB) + switch (ktype) { + case TYPE_SB: typeData[3] = A4_SB; + break; + case TYPE_OOP: + sortOop(typeData); + break; + } key = MinObject.getKey(type, typeData[0], typeData[1], typeData[2], typeData[3]); if (ffParams.containsKey(key)) @@ -838,20 +847,14 @@ return key; break; case TYPE_TORSION: - key = MinObject.getKey(type, getEquivalentType(typeData[0], 0), - typeData[1], typeData[2], getEquivalentType(typeData[3], 2)); - if (ffParams.containsKey(key)) + if (ffParams.containsKey(key = getTorsionKey(type, 0, 2))) return key; - key = MinObject.getKey(type, getEquivalentType(typeData[0], 2), - typeData[1], typeData[2], getEquivalentType(typeData[3], 0)); - if (ffParams.containsKey(key)) + if (ffParams.containsKey(key = getTorsionKey(type, 2, 0))) return key; - key = MinObject.getKey(type, getEquivalentType(typeData[0], 2), - typeData[1], typeData[2], getEquivalentType(typeData[3], 2)); + if (ffParams.containsKey(key = getTorsionKey(type, 2, 2))) + return key; // set type = 0 if necessary... - return (ffParams.containsKey(key) ? key : MinObject.getKey(0, - getEquivalentType(typeData[0], 2), typeData[1], typeData[2], - getEquivalentType(typeData[3], 2))); + return getTorsionKey(0, 2, 2); case TYPE_SB: // use periodic row info int r1 = getRowFor(data[0]); @@ -873,6 +876,9 @@ case TYPE_ANGLE: isSwapped = (fixTypeOrder(typeData, 0, 2)); break; + case TYPE_OOP: + sortOop(typeData); + break; } key = MinObject.getKey(type, typeData[0], typeData[1], typeData[2], typeData[3]); @@ -893,6 +899,11 @@ return key; } + private Integer getTorsionKey(int type, int i, int j) { + return MinObject.getKey(type, getEquivalentType(typeData[0], i), + typeData[1], typeData[2], getEquivalentType(typeData[3], j)); + } + private Integer applyEmpiricalRules(MinObject o, int ktype) { // TODO return null; @@ -906,29 +917,15 @@ private int[] typeData = new int[4]; double getOutOfPlaneParameter(int[] data) { - int typeA = typeData[0] = typeOf(data[0]); - /*int typeB = */typeData[1] = typeOf(data[1]); - int typeC = typeData[2] = typeOf(data[2]); - int typeD = typeData[3] = typeOf(data[3]); - Object params = getSortedOOPParam(typeData); - if (params == null) { - for (int i = 0; i < 3; i++) { - typeData[0] = getEquivalentType(typeA, i); - typeData[2] = getEquivalentType(typeC, i); - typeData[3] = getEquivalentType(typeD, i); - if ((params = getSortedOOPParam(typeData)) != null) - break; - } - } - return ((Double) params).doubleValue(); + Integer key = getKey(data, KEY_OOP, TYPE_OOP); + Double params = (Double) ffParams.get(key); + return (params == null ? 0 : params.doubleValue()); } - private static Object getSortedOOPParam(int[] typeData) { + private static void sortOop(int[] typeData) { fixTypeOrder(typeData, 0, 2); fixTypeOrder(typeData, 0, 3); fixTypeOrder(typeData, 2, 3); - return ffParams.get(MinObject.getKey(KEY_OOP, - typeData[0], typeData[1], typeData[2], typeData[3])); } /** Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/mmff/validate/checkmm.spt =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/mmff/validate/checkmm.spt 2012-05-17 15:24:00 UTC (rev 17151) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/mmff/validate/checkmm.spt 2012-05-17 16:10:31 UTC (rev 17152) @@ -1,3 +1,4 @@ +set forcefield "MMFF" print "checkmm.spt reading mmff-validation.dat..." vdat = load("mmff-validation.dat").split("model ") tdat = {} @@ -160,6 +161,9 @@ set useMinimizationThread false set minimizationSilent var nBad = 0; + var a = []; + var max = 0.1; + a[last + 1] = 0; for (var m = 0; m <= last; m++) { model @{m+1} select {modelindex=m} @@ -167,10 +171,12 @@ var x = 0.0 + getproperty("minimizationinfo").split("TOTAL ENERGY =")[2].trim().split(" ")[1] var ref = mmdat[_modelName]["energy"] var diff = abs(x - ref) - if (diff > 0.1) { + a[m + 1] = diff + if (diff > max) { print "" + (++nBad) + "\t" + _modelName + "\t E= \t" + x + "\t Eref=\t" + ref + "\t diff= \t" + diff } } + print "for " + (last+1) + " atoms, " + nBad + " have energies differences outside the range " + (-max) + " to " + max + " with a standard deviation of " + a.stddev } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2014-02-25 00:33:33
|
Revision: 19381 http://sourceforge.net/p/jmol/code/19381 Author: hansonr Date: 2014-02-25 00:33:13 +0000 (Tue, 25 Feb 2014) Log Message: ----------- stretch-bend report error Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/MinObject.java trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java Modified: trunk/Jmol/src/org/jmol/minimize/MinObject.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/MinObject.java 2014-02-24 22:51:20 UTC (rev 19380) +++ trunk/Jmol/src/org/jmol/minimize/MinObject.java 2014-02-25 00:33:13 UTC (rev 19381) @@ -33,11 +33,11 @@ int c = i & 0x7F; i >>= 7; int d = i & 0x7F; - return type + ": " + return (type < 0 ? type + ": " : "") + (a < 10 ? " " : " ") + a + (b < 10 ? " " : " ") + b + (c < 10 ? " " : " ") + c - + (d < 10 ? " " : " ") + d; + + (d > 120 ? "" : (d < 10 ? " " : " ") + d); } } Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java 2014-02-24 22:51:20 UTC (rev 19380) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/CalculationsMMFF.java 2014-02-25 00:33:13 UTC (rev 19381) @@ -306,13 +306,16 @@ double theta0 = datakat0[1]; double r0ij = dataij[1]; double r0jk = datajk[1]; - calc.addLast(new Object[] { angle.data, new double[] { data[0], theta0, r0ij } }); + calc.addLast(new Object[] { angle.data, new double[] { data[0], theta0, r0ij }, angle.sbKey }); calc.addLast(new Object[] { new int[] {angle.data[2], angle.data[1], angle.data[0]}, - new double[] { data[1], theta0, r0jk } }); + new double[] { data[1], theta0, r0jk }, angle.sbKey }); } @Override double compute(Object[] dataIn) { + + key = (Integer) dataIn[2]; + getPointers(dataIn); double k = 2.51210 * dData[0]; double t0 = dData[1]; @@ -589,7 +592,7 @@ case CALC_ANGLE: case CALC_STRETCH_BEND: return Txt.sprintf( - "%15s %-5s %-5s %-5s %8.3f %8.3f %8.3f %8.3f", + "%11s %-5s %-5s %-5s %8.3f %8.3f %8.3f %8.3f", "ssssFI", new Object[] { MinObject.decodeKey(c.key), minAtoms[c.ia].sType, minAtoms[c.ib].sType, minAtoms[c.ic].sType, new float[] { (float)(c.theta * RAD_TO_DEG), (float) c.dData[1] /*THETA0*/, This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <ha...@us...> - 2023-07-02 22:44:21
|
Revision: 22525 http://sourceforge.net/p/jmol/code/22525 Author: hansonr Date: 2023-07-02 22:44:19 +0000 (Sun, 02 Jul 2023) Log Message: ----------- Jmol.___JmolVersion="16.1.17" // (legacy) also 16.1.18 (swingJS) bug fix: set picking assignatom_P (specifically phosphorus) fails (second try!) bug fix: MINIMIZE crystal structure better optimization bug fix: LOAD APPEND VAR does not work bug fix: LOAD "" should return last actual file loaded, not "string" if LOAD VAR or LOAD INLINE was used new feature: MODELKIT MINIMIZE - for crystal structures in the current model only - does not need anything more than the asymmetric unit atoms - carries out a minimization on the asymmetric unit within a 27-cell block {444 666 1} - applies that minimization to all atoms in the current model - temporarily creates a new CIF model, minimizes it, then deletes it - can accept a limited number of additional parameters such as SILENT and - Equivalent to function minimizeXtal() { if (within("unitCell") == 0) return; var m = {thisModel}; var au = {thisModel & !symmetry}; var x = write("cif"); var usethread = useMinimizationThread; var apnew = appendNew; set useMinimizationThread false; set appendNew true; load append var x {444 666 1}; frame last; minimize; set useMinimizationThread @usethread; set appendNew @apnew; var pts = {thisModel & symop=1555}.xyz.all; zap {thisModel}; frame @m; modelkit moveto @au @pts; } new feature: MODELKIT MOVETO @atoms @points - similar to MODELKIT MOVETO @atom @xyz, but for a set of atoms - atoms may be just the asymmetric unit - points is a set of points, same length as atoms - bulk conversion of a set of atom positions, applying symmetry Modified Paths: -------------- trunk/Jmol/src/org/jmol/minimize/Minimizer.java trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java Modified: trunk/Jmol/src/org/jmol/minimize/Minimizer.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2023-07-02 20:53:18 UTC (rev 22524) +++ trunk/Jmol/src/org/jmol/minimize/Minimizer.java 2023-07-02 22:44:19 UTC (rev 22525) @@ -169,7 +169,7 @@ try { setEnergyUnits(); - // if the user indicated minimize ... FIX ... or we don't have any defualt, + // if the user indicated minimize ... FIX ... or we don't have any default, // use the bsFixed coming in here, which is set to "nearby and in frame" in that case. // and if something is fixed, then AND it with "nearby and in frame" as well. if (!haveFixed && bsFixedDefault != null) Modified: trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java =================================================================== --- trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java 2023-07-02 20:53:18 UTC (rev 22524) +++ trunk/Jmol/src/org/jmol/minimize/forcefield/ForceField.java 2023-07-02 22:44:19 UTC (rev 22525) @@ -209,7 +209,6 @@ } return false; } - //System.out.println(isPreliminary + " " + getNormalizedDE() + " " + currentStep); if (isPreliminary && getNormalizedDE() >= 2) { // looking back at this after some time, I don't exactly see why I wanted // this to stay in preliminary mode unless |DE| >= 2 * crit. @@ -262,8 +261,6 @@ atom.force[0] = -getDE(atom, terms, 0, delta); atom.force[1] = -getDE(atom, terms, 1, delta); atom.force[2] = -getDE(atom, terms, 2, delta); - //if (atom.atom.getAtomIndex() == 2) - //System.out.println(" atom + " + atom.atom.getAtomIndex() + " force=" + atom.force[0] + " " + atom.force[1] + " " + atom.force[2] ); return; } @@ -387,10 +384,10 @@ double step = 0.75 * trustRadius; double trustRadius2 = trustRadius * trustRadius; - double e1 = energyFull(false, true); int nSteps = 10; - for (int iStep = 0; iStep < nSteps; iStep++) { + boolean isDone = false; + for (int iStep = 0; iStep < nSteps && !isDone; iStep++) { saveCoordinates(); for (int i = 0; i < minAtomCount; ++i) { if (bsMinFixed == null || !bsMinFixed.get(i)) { @@ -398,8 +395,9 @@ double[] coord = minAtoms[i].coord; double f2 = (force[0] * force[0] + force[1] * force[1] + force[2] * force[2]); - if (f2 > trustRadius2 / step / step) { - f2 = trustRadius / Math.sqrt(f2) / step; + double f = trustRadius2 / step / step / f2; + if (1 > f) { + f2 = Math.sqrt(f); force[0] *= f2; force[1] *= f2; force[2] *= f2; @@ -417,11 +415,12 @@ } } } - if (doUpdateAtoms) + if (doUpdateAtoms) { + // only necessary if symmetry constraints might adjust positions minimizer.updateAtomXYZ(); + } double e2 = energyFull(false, true); - if (Util.isNear3(e2, e1, 1.0e-3)) - break; + isDone = Util.isNear3(e2, e1, 1.0e-3); if (e2 > e1) { step *= 0.1; restoreCoordinates(); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |