diff --git a/examples/QMMM/.gitkeep b/examples/QMMM/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/examples/QMMM/2E4E.pdb b/examples/QMMM/2E4E.pdb new file mode 100644 index 0000000..ac2283c --- /dev/null +++ b/examples/QMMM/2E4E.pdb @@ -0,0 +1,292 @@ +HEADER DE NOVO PROTEIN 06-DEC-06 2E4E +TITLE NMR STRUCTURE OF D4P/K7G MUTANT OF GPM12 +COMPND MOL_ID: 1; +COMPND 2 MOLECULE: GPM12; +COMPND 3 CHAIN: A; +COMPND 4 SYNONYM: CHIGNOLIN; +COMPND 5 ENGINEERED: YES; +COMPND 6 MUTATION: YES +SOURCE MOL_ID: 1; +SOURCE 2 SYNTHETIC: YES; +SOURCE 3 OTHER_DETAILS: CHEMICAL PEPTIDE SYNTHESIS +KEYWDS BETA-HAIRPIN, MINI-PROTEIN, CHIGNOLIN, B1 DOMAIN OF PROTEIN G, DE +KEYWDS 2 NOVO PROTEIN +EXPDTA SOLUTION NMR +NUMMDL 23 +AUTHOR T.TERADA,D.SATOH,T.MIKAWA,Y.ITO,K.SHIMIZU +REVDAT 4 29-MAY-24 2E4E 1 REMARK +REVDAT 3 09-MAR-22 2E4E 1 REMARK +REVDAT 2 24-FEB-09 2E4E 1 VERSN +REVDAT 1 05-FEB-08 2E4E 0 +JRNL AUTH T.TERADA,D.SATOH,T.MIKAWA,Y.ITO,K.SHIMIZU +JRNL TITL UNDERSTANDING THE ROLES OF AMINO ACID RESIDUES IN TERTIARY +JRNL TITL 2 STRUCTURE FORMATION OF CHIGNOLIN BY USING MOLECULAR DYNAMICS +JRNL TITL 3 SIMULATION +JRNL REF TO BE PUBLISHED +JRNL REFN +REMARK 2 +REMARK 2 RESOLUTION. NOT APPLICABLE. +REMARK 3 +REMARK 3 REFINEMENT. +REMARK 3 PROGRAM : XWINNMR 3.5, CNS 1.1 +REMARK 3 AUTHORS : BRUKER BIOSPIN CORPORATION (XWINNMR), +REMARK 3 BRUNGER,ADAMS,CLORE,DELANO,GROS,GROSSE-KUNSTLEVE,JIANG,KUSZEWSKI, +REMARK 3 NILGES, PANNU,READ,RICE,SIMONSON,WARREN (CNS) +REMARK 3 +REMARK 3 OTHER REFINEMENT REMARKS: THE STRUCTURES ARE BASED ON 119 NOE +REMARK 3 -DERIVED DISTANCE CONSTRAINTS. +REMARK 4 +REMARK 4 2E4E COMPLIES WITH FORMAT V. 3.15, 01-DEC-08 +REMARK 100 +REMARK 100 THIS ENTRY HAS BEEN PROCESSED BY PDBJ ON 13-DEC-06. +REMARK 100 THE DEPOSITION ID IS D_1000026210. +REMARK 210 +REMARK 210 EXPERIMENTAL DETAILS +REMARK 210 EXPERIMENT TYPE : NMR +REMARK 210 TEMPERATURE (KELVIN) : 277 +REMARK 210 PH : 5.5 +REMARK 210 IONIC STRENGTH : NULL +REMARK 210 PRESSURE : 1 ATM +REMARK 210 SAMPLE CONTENTS : 2MM GPM12(D4P/K7G) +REMARK 210 +REMARK 210 NMR EXPERIMENTS CONDUCTED : 2D NOESY; 2D TOCSY; 2D ROESY +REMARK 210 SPECTROMETER FIELD STRENGTH : 600 MHZ +REMARK 210 SPECTROMETER MODEL : DRX +REMARK 210 SPECTROMETER MANUFACTURER : BRUKER +REMARK 210 +REMARK 210 STRUCTURE DETERMINATION. +REMARK 210 SOFTWARE USED : AZARA 2.7, ANSIG 3.3 FOR OPENGL +REMARK 210 VERSION 1.0.6, CNS 1.1 +REMARK 210 METHOD USED : SIMULATED ANNEALING +REMARK 210 +REMARK 210 CONFORMERS, NUMBER CALCULATED : 200 +REMARK 210 CONFORMERS, NUMBER SUBMITTED : 23 +REMARK 210 CONFORMERS, SELECTION CRITERIA : STRUCTURES WITH THE LEAST +REMARK 210 RESTRAINT VIOLATIONS +REMARK 210 +REMARK 210 BEST REPRESENTATIVE CONFORMER IN THIS ENSEMBLE : 1 +REMARK 210 +REMARK 210 REMARK: THIS STRUCTURE WAS DETERMINED USING STANDARD 2D +REMARK 210 HOMONUCLEAR TECHNIQUES. +REMARK 215 +REMARK 215 NMR STUDY +REMARK 215 THE COORDINATES IN THIS ENTRY WERE GENERATED FROM SOLUTION +REMARK 215 NMR DATA. PROTEIN DATA BANK CONVENTIONS REQUIRE THAT +REMARK 215 CRYST1 AND SCALE RECORDS BE INCLUDED, BUT THE VALUES ON +REMARK 215 THESE RECORDS ARE MEANINGLESS. +REMARK 300 +REMARK 300 BIOMOLECULE: 1 +REMARK 300 SEE REMARK 350 FOR THE AUTHOR PROVIDED AND/OR PROGRAM +REMARK 300 GENERATED ASSEMBLY INFORMATION FOR THE STRUCTURE IN +REMARK 300 THIS ENTRY. THE REMARK MAY ALSO PROVIDE INFORMATION ON +REMARK 300 BURIED SURFACE AREA. +REMARK 350 +REMARK 350 COORDINATES FOR A COMPLETE MULTIMER REPRESENTING THE KNOWN +REMARK 350 BIOLOGICALLY SIGNIFICANT OLIGOMERIZATION STATE OF THE +REMARK 350 MOLECULE CAN BE GENERATED BY APPLYING BIOMT TRANSFORMATIONS +REMARK 350 GIVEN BELOW. BOTH NON-CRYSTALLOGRAPHIC AND +REMARK 350 CRYSTALLOGRAPHIC OPERATIONS ARE GIVEN. +REMARK 350 +REMARK 350 BIOMOLECULE: 1 +REMARK 350 AUTHOR DETERMINED BIOLOGICAL UNIT: MONOMERIC +REMARK 350 APPLY THE FOLLOWING TO CHAINS: A +REMARK 350 BIOMT1 1 1.000000 0.000000 0.000000 0.00000 +REMARK 350 BIOMT2 1 0.000000 1.000000 0.000000 0.00000 +REMARK 350 BIOMT3 1 0.000000 0.000000 1.000000 0.00000 +REMARK 500 +REMARK 500 GEOMETRY AND STEREOCHEMISTRY +REMARK 500 SUBTOPIC: TORSION ANGLES +REMARK 500 +REMARK 500 TORSION ANGLES OUTSIDE THE EXPECTED RAMACHANDRAN REGIONS: +REMARK 500 (M=MODEL NUMBER; RES=RESIDUE NAME; C=CHAIN IDENTIFIER; +REMARK 500 SSEQ=SEQUENCE NUMBER; I=INSERTION CODE). +REMARK 500 +REMARK 500 STANDARD TABLE: +REMARK 500 FORMAT:(10X,I3,1X,A3,1X,A1,I4,A1,4X,F7.2,3X,F7.2) +REMARK 500 +REMARK 500 EXPECTED VALUES: GJ KLEYWEGT AND TA JONES (1996). PHI/PSI- +REMARK 500 CHOLOGY: RAMACHANDRAN REVISITED. STRUCTURE 4, 1395 - 1400 +REMARK 500 +REMARK 500 M RES CSSEQI PSI PHI +REMARK 500 1 THR A 8 -137.71 -113.44 +REMARK 500 1 PHE A 9 10.41 -141.68 +REMARK 500 2 THR A 8 -138.60 -119.39 +REMARK 500 3 THR A 8 -138.18 -119.00 +REMARK 500 4 THR A 8 -137.16 -115.24 +REMARK 500 5 THR A 8 -137.68 -115.90 +REMARK 500 5 PHE A 9 10.11 -142.20 +REMARK 500 6 THR A 8 -138.24 -112.61 +REMARK 500 7 THR A 8 -138.24 -113.21 +REMARK 500 7 PHE A 9 10.49 -140.88 +REMARK 500 8 THR A 8 -137.39 -119.00 +REMARK 500 9 THR A 8 -137.53 -116.31 +REMARK 500 9 PHE A 9 10.58 -140.69 +REMARK 500 10 THR A 8 -137.61 -119.71 +REMARK 500 11 THR A 8 -137.28 -121.13 +REMARK 500 12 THR A 8 -137.20 -118.51 +REMARK 500 13 THR A 8 -137.46 -118.27 +REMARK 500 14 THR A 8 -137.55 -123.85 +REMARK 500 15 THR A 8 -138.05 -116.17 +REMARK 500 15 PHE A 9 10.46 -140.38 +REMARK 500 16 THR A 8 -137.49 -117.34 +REMARK 500 17 THR A 8 -138.57 -116.29 +REMARK 500 18 THR A 8 -137.28 -117.18 +REMARK 500 18 PHE A 9 10.01 -141.14 +REMARK 500 19 THR A 8 -138.00 -113.18 +REMARK 500 20 THR A 8 -137.18 -119.00 +REMARK 500 21 THR A 8 -137.86 -110.95 +REMARK 500 21 PHE A 9 10.25 -144.25 +REMARK 500 22 THR A 8 -137.90 -115.13 +REMARK 500 22 PHE A 9 10.23 -140.13 +REMARK 500 23 THR A 8 -137.53 -120.32 +REMARK 500 +REMARK 500 REMARK: NULL +REMARK 999 +REMARK 999 SEQUENCE +REMARK 999 THERE ARE MUTANTS D4P, K7G. THESE MUTATIONS CONVERT +REMARK 999 THE DISORDERD STRUCTURE OF GPM12 INTO A CHIGNOLIN-LIKE +REMARK 999 ORDERED STRUCTURE. +REMARK 999 SEQUENCE OF RESIDUES 2-9 OF GPM12 IS THE SAME AS +REMARK 999 THAT OF RESIDUES 45-52 OF THE B1 DOMAIN OF PROTEIN G. +DBREF 2E4E A 1 10 PDB 2E4E 2E4E 1 10 +SEQRES 1 A 10 GLY TYR ASP PRO ALA THR GLY THR PHE GLY +CRYST1 1.000 1.000 1.000 90.00 90.00 90.00 P 1 1 +ORIGX1 1.000000 0.000000 0.000000 0.00000 +ORIGX2 0.000000 1.000000 0.000000 0.00000 +ORIGX3 0.000000 0.000000 1.000000 0.00000 +SCALE1 1.000000 0.000000 0.000000 0.00000 +SCALE2 0.000000 1.000000 0.000000 0.00000 +SCALE3 0.000000 0.000000 1.000000 0.00000 +MODEL 1 +ATOM 1 N GLY A 1 -0.381 -7.452 0.697 1.00 0.00 N +ATOM 2 CA GLY A 1 -1.570 -6.645 0.308 1.00 0.00 C +ATOM 3 C GLY A 1 -1.307 -5.155 0.383 1.00 0.00 C +ATOM 4 O GLY A 1 -1.508 -4.433 -0.594 1.00 0.00 O +ATOM 5 H1 GLY A 1 0.285 -7.516 -0.100 1.00 0.00 H +ATOM 6 H2 GLY A 1 -0.672 -8.413 0.968 1.00 0.00 H +ATOM 7 H3 GLY A 1 0.105 -7.010 1.504 1.00 0.00 H +ATOM 8 HA2 GLY A 1 -1.849 -6.900 -0.703 1.00 0.00 H +ATOM 9 HA3 GLY A 1 -2.388 -6.890 0.969 1.00 0.00 H +ATOM 10 N TYR A 2 -0.857 -4.692 1.545 1.00 0.00 N +ATOM 11 CA TYR A 2 -0.566 -3.279 1.743 1.00 0.00 C +ATOM 12 C TYR A 2 0.530 -2.813 0.797 1.00 0.00 C +ATOM 13 O TYR A 2 1.589 -3.430 0.707 1.00 0.00 O +ATOM 14 CB TYR A 2 -0.139 -3.019 3.189 1.00 0.00 C +ATOM 15 CG TYR A 2 0.174 -1.565 3.473 1.00 0.00 C +ATOM 16 CD1 TYR A 2 -0.841 -0.667 3.771 1.00 0.00 C +ATOM 17 CD2 TYR A 2 1.481 -1.090 3.437 1.00 0.00 C +ATOM 18 CE1 TYR A 2 -0.566 0.663 4.026 1.00 0.00 C +ATOM 19 CE2 TYR A 2 1.765 0.236 3.691 1.00 0.00 C +ATOM 20 CZ TYR A 2 0.740 1.109 3.985 1.00 0.00 C +ATOM 21 OH TYR A 2 1.019 2.433 4.238 1.00 0.00 O +ATOM 22 H TYR A 2 -0.717 -5.314 2.285 1.00 0.00 H +ATOM 23 HA TYR A 2 -1.467 -2.720 1.538 1.00 0.00 H +ATOM 24 HB2 TYR A 2 -0.937 -3.322 3.853 1.00 0.00 H +ATOM 25 HB3 TYR A 2 0.744 -3.599 3.407 1.00 0.00 H +ATOM 26 HD1 TYR A 2 -1.860 -1.019 3.800 1.00 0.00 H +ATOM 27 HD2 TYR A 2 2.284 -1.773 3.207 1.00 0.00 H +ATOM 28 HE1 TYR A 2 -1.369 1.346 4.257 1.00 0.00 H +ATOM 29 HE2 TYR A 2 2.788 0.583 3.656 1.00 0.00 H +ATOM 30 HH TYR A 2 1.827 2.498 4.753 1.00 0.00 H +ATOM 31 N ASP A 3 0.268 -1.712 0.103 1.00 0.00 N +ATOM 32 CA ASP A 3 1.229 -1.141 -0.828 1.00 0.00 C +ATOM 33 C ASP A 3 2.107 -0.122 -0.107 1.00 0.00 C +ATOM 34 O ASP A 3 1.659 0.983 0.196 1.00 0.00 O +ATOM 35 CB ASP A 3 0.505 -0.473 -1.998 1.00 0.00 C +ATOM 36 CG ASP A 3 1.338 -0.464 -3.264 1.00 0.00 C +ATOM 37 OD1 ASP A 3 2.286 -1.271 -3.360 1.00 0.00 O +ATOM 38 OD2 ASP A 3 1.042 0.353 -4.164 1.00 0.00 O +ATOM 39 H ASP A 3 -0.588 -1.269 0.228 1.00 0.00 H +ATOM 40 HA ASP A 3 1.840 -1.941 -1.200 1.00 0.00 H +ATOM 41 HB2 ASP A 3 -0.413 -1.008 -2.198 1.00 0.00 H +ATOM 42 HB3 ASP A 3 0.272 0.548 -1.733 1.00 0.00 H +ATOM 43 N PRO A 4 3.363 -0.485 0.210 1.00 0.00 N +ATOM 44 CA PRO A 4 4.283 0.397 0.932 1.00 0.00 C +ATOM 45 C PRO A 4 4.825 1.530 0.068 1.00 0.00 C +ATOM 46 O PRO A 4 5.087 2.629 0.561 1.00 0.00 O +ATOM 47 CB PRO A 4 5.424 -0.534 1.362 1.00 0.00 C +ATOM 48 CG PRO A 4 4.908 -1.903 1.123 1.00 0.00 C +ATOM 49 CD PRO A 4 3.992 -1.776 -0.066 1.00 0.00 C +ATOM 50 HA PRO A 4 3.815 0.802 1.808 1.00 0.00 H +ATOM 51 HB2 PRO A 4 6.303 -0.334 0.771 1.00 0.00 H +ATOM 52 HB3 PRO A 4 5.635 -0.376 2.409 1.00 0.00 H +ATOM 53 HG2 PRO A 4 5.727 -2.584 0.929 1.00 0.00 H +ATOM 54 HG3 PRO A 4 4.351 -2.218 1.988 1.00 0.00 H +ATOM 55 HD2 PRO A 4 4.548 -1.746 -1.000 1.00 0.00 H +ATOM 56 HD3 PRO A 4 3.267 -2.569 -0.074 1.00 0.00 H +ATOM 57 N ALA A 5 4.995 1.258 -1.220 1.00 0.00 N +ATOM 58 CA ALA A 5 5.510 2.257 -2.149 1.00 0.00 C +ATOM 59 C ALA A 5 4.603 3.476 -2.189 1.00 0.00 C +ATOM 60 O ALA A 5 5.046 4.603 -1.961 1.00 0.00 O +ATOM 61 CB ALA A 5 5.657 1.660 -3.540 1.00 0.00 C +ATOM 62 H ALA A 5 4.769 0.368 -1.552 1.00 0.00 H +ATOM 63 HA ALA A 5 6.490 2.558 -1.805 1.00 0.00 H +ATOM 64 HB1 ALA A 5 4.776 1.881 -4.121 1.00 0.00 H +ATOM 65 HB2 ALA A 5 5.779 0.590 -3.461 1.00 0.00 H +ATOM 66 HB3 ALA A 5 6.524 2.087 -4.024 1.00 0.00 H +ATOM 67 N THR A 6 3.330 3.242 -2.476 1.00 0.00 N +ATOM 68 CA THR A 6 2.353 4.302 -2.546 1.00 0.00 C +ATOM 69 C THR A 6 1.744 4.558 -1.167 1.00 0.00 C +ATOM 70 O THR A 6 1.312 5.669 -0.865 1.00 0.00 O +ATOM 71 CB THR A 6 1.272 3.924 -3.556 1.00 0.00 C +ATOM 72 OG1 THR A 6 1.778 3.992 -4.878 1.00 0.00 O +ATOM 73 CG2 THR A 6 0.053 4.805 -3.493 1.00 0.00 C +ATOM 74 H THR A 6 3.037 2.330 -2.646 1.00 0.00 H +ATOM 75 HA THR A 6 2.857 5.191 -2.881 1.00 0.00 H +ATOM 76 HB THR A 6 0.957 2.908 -3.365 1.00 0.00 H +ATOM 77 HG1 THR A 6 1.480 3.224 -5.374 1.00 0.00 H +ATOM 78 HG21 THR A 6 -0.252 4.915 -2.465 1.00 0.00 H +ATOM 79 HG22 THR A 6 -0.746 4.352 -4.058 1.00 0.00 H +ATOM 80 HG23 THR A 6 0.288 5.773 -3.907 1.00 0.00 H +ATOM 81 N GLY A 7 1.720 3.519 -0.337 1.00 0.00 N +ATOM 82 CA GLY A 7 1.174 3.640 0.995 1.00 0.00 C +ATOM 83 C GLY A 7 -0.327 3.454 1.037 1.00 0.00 C +ATOM 84 O GLY A 7 -1.019 4.082 1.839 1.00 0.00 O +ATOM 85 H GLY A 7 2.083 2.662 -0.630 1.00 0.00 H +ATOM 86 HA2 GLY A 7 1.626 2.889 1.620 1.00 0.00 H +ATOM 87 HA3 GLY A 7 1.419 4.613 1.380 1.00 0.00 H +ATOM 88 N THR A 8 -0.827 2.582 0.175 1.00 0.00 N +ATOM 89 CA THR A 8 -2.212 2.286 0.091 1.00 0.00 C +ATOM 90 C THR A 8 -2.445 0.850 0.525 1.00 0.00 C +ATOM 91 O THR A 8 -1.870 0.377 1.503 1.00 0.00 O +ATOM 92 CB THR A 8 -2.665 2.525 -1.349 1.00 0.00 C +ATOM 93 OG1 THR A 8 -1.848 3.498 -1.976 1.00 0.00 O +ATOM 94 CG2 THR A 8 -4.084 3.001 -1.440 1.00 0.00 C +ATOM 95 H THR A 8 -0.247 2.117 -0.430 1.00 0.00 H +ATOM 96 HA THR A 8 -2.749 2.936 0.740 1.00 0.00 H +ATOM 97 HB THR A 8 -2.582 1.605 -1.906 1.00 0.00 H +ATOM 98 HG1 THR A 8 -1.921 4.331 -1.503 1.00 0.00 H +ATOM 99 HG21 THR A 8 -4.091 4.057 -1.656 1.00 0.00 H +ATOM 100 HG22 THR A 8 -4.574 2.822 -0.499 1.00 0.00 H +ATOM 101 HG23 THR A 8 -4.591 2.464 -2.224 1.00 0.00 H +ATOM 102 N PHE A 9 -3.288 0.178 -0.207 1.00 0.00 N +ATOM 103 CA PHE A 9 -3.628 -1.203 0.064 1.00 0.00 C +ATOM 104 C PHE A 9 -3.786 -2.000 -1.216 1.00 0.00 C +ATOM 105 O PHE A 9 -4.256 -3.138 -1.216 1.00 0.00 O +ATOM 106 CB PHE A 9 -4.913 -1.242 0.854 1.00 0.00 C +ATOM 107 CG PHE A 9 -4.708 -1.193 2.342 1.00 0.00 C +ATOM 108 CD1 PHE A 9 -4.643 0.023 3.004 1.00 0.00 C +ATOM 109 CD2 PHE A 9 -4.577 -2.360 3.076 1.00 0.00 C +ATOM 110 CE1 PHE A 9 -4.452 0.073 4.372 1.00 0.00 C +ATOM 111 CE2 PHE A 9 -4.385 -2.316 4.443 1.00 0.00 C +ATOM 112 CZ PHE A 9 -4.322 -1.099 5.091 1.00 0.00 C +ATOM 113 H PHE A 9 -3.695 0.626 -0.946 1.00 0.00 H +ATOM 114 HA PHE A 9 -2.836 -1.621 0.631 1.00 0.00 H +ATOM 115 HB2 PHE A 9 -5.504 -0.392 0.561 1.00 0.00 H +ATOM 116 HB3 PHE A 9 -5.444 -2.142 0.612 1.00 0.00 H +ATOM 117 HD1 PHE A 9 -4.744 0.939 2.442 1.00 0.00 H +ATOM 118 HD2 PHE A 9 -4.626 -3.312 2.568 1.00 0.00 H +ATOM 119 HE1 PHE A 9 -4.403 1.026 4.877 1.00 0.00 H +ATOM 120 HE2 PHE A 9 -4.283 -3.233 5.003 1.00 0.00 H +ATOM 121 HZ PHE A 9 -4.173 -1.064 6.161 1.00 0.00 H +ATOM 122 N GLY A 10 -3.386 -1.377 -2.293 1.00 0.00 N +ATOM 123 CA GLY A 10 -3.467 -1.996 -3.604 1.00 0.00 C +ATOM 124 C GLY A 10 -3.277 -0.996 -4.728 1.00 0.00 C +ATOM 125 O GLY A 10 -2.140 -0.507 -4.901 1.00 0.00 O +ATOM 126 OXT GLY A 10 -4.264 -0.703 -5.435 1.00 0.00 O +ATOM 127 H GLY A 10 -3.031 -0.480 -2.191 1.00 0.00 H +ATOM 128 HA2 GLY A 10 -2.703 -2.755 -3.678 1.00 0.00 H +ATOM 129 HA3 GLY A 10 -4.436 -2.462 -3.711 1.00 0.00 H +TER 130 GLY A 10 +ENDMDL +END diff --git a/examples/QMMM/2E4E_RHF-DFT-QMMM_energy.inp b/examples/QMMM/2E4E_RHF-DFT-QMMM_energy.inp new file mode 100644 index 0000000..919340b --- /dev/null +++ b/examples/QMMM/2E4E_RHF-DFT-QMMM_energy.inp @@ -0,0 +1,19 @@ +[input] +system=2E4E.pdb 1 2 5 6 7 8 9 3 4 10 22 11 23 +charge=1 +runtype=energy +functional=bhhlyp +basis=6-31gs +method=hf +qmmm_flag=True + +[guess] +type=huckel + +[scf] +multiplicity=1 +type=rhf + +[dftgrid] +rad_type=becke + diff --git a/examples/QMMM/ala.inp b/examples/QMMM/ala.inp new file mode 100644 index 0000000..9a30dd7 --- /dev/null +++ b/examples/QMMM/ala.inp @@ -0,0 +1,19 @@ +[input] +system=ala.pdb 9 10 17 18 19 +charge=0 +runtype=energy +functional=bhhlyp +basis=6-31gs +method=hf +qmmm_flag=True + +[guess] +type=huckel + +[scf] +multiplicity=1 +type=rhf + +[dftgrid] +rad_type=becke + diff --git a/examples/QMMM/ala.pdb b/examples/QMMM/ala.pdb new file mode 100644 index 0000000..c723fc5 --- /dev/null +++ b/examples/QMMM/ala.pdb @@ -0,0 +1,23 @@ +HEADER ala +COMPND +SOURCE +ATOM 1 CH3 ACE 1 -0.028 0.101 -0.091 +ATOM 2 C ACE 1 -0.023 0.157 1.419 +ATOM 3 O ACE 1 1.014 0.404 2.021 +ATOM 4 H1 ACE 1 0.978 0.305 -0.455 +ATOM 5 H2 ACE 1 -0.715 0.851 -0.477 +ATOM 6 H3 ACE 1 -0.339 -0.891 -0.414 +ATOM 7 N ALA 2 -1.183 -0.077 2.026 +ATOM 8 CA ALA 2 -1.370 -0.108 3.477 +ATOM 9 C ALA 2 -2.450 -1.132 3.872 +ATOM 10 O ALA 2 -3.296 -1.508 3.068 +ATOM 11 H ALA 2 -1.986 -0.315 1.463 +ATOM 12 HA ALA 2 -0.433 -0.408 3.951 +ATOM 13 CB ALA 2 -1.730 1.305 3.959 +ATOM 14 HB1 ALA 2 -1.842 1.311 5.044 +ATOM 15 HB2 ALA 2 -2.665 1.629 3.500 +ATOM 16 HB3 ALA 2 -0.935 2.000 3.685 +ATOM 17 N NH2 3 -2.452 -1.588 5.119 +ATOM 18 H1 NH2 3 -1.760 -1.272 5.777 +ATOM 19 H2 NH2 3 -3.162 -2.253 5.377 +END diff --git a/include/oqp.h b/include/oqp.h index 5024dfa..14f90b3 100644 --- a/include/oqp.h +++ b/include/oqp.h @@ -77,7 +77,7 @@ struct dft_parameters { bool dft_wt_der; }; - struct tddft_parameters { +struct tddft_parameters { int64_t nstate; int64_t target_state; int64_t maxvec; @@ -123,6 +123,7 @@ struct control_parameters { double esp_constr; bool basis_set_issue; double conf_print_threshold; + bool qmmm_flag; }; struct mpi_communicator { diff --git a/pyoqp/oqp/library/ints_1e.py b/pyoqp/oqp/library/ints_1e.py index 8644f13..04bb260 100644 --- a/pyoqp/oqp/library/ints_1e.py +++ b/pyoqp/oqp/library/ints_1e.py @@ -1,5 +1,11 @@ import oqp +import numpy as np +from oqp.utils.qmmm import openmm_potential,openmm_energy def ints_1e(mol): """Compute a set of one-electron integrals""" + if mol.config['input']['qmmm_flag']: + current_xyz = mol.get_system().reshape((-1, 3)) + mol.data["OQP::mm_potential"]=np.array(openmm_potential(current_xyz)).tolist() + mol.data["OQP::mm_energy"]=openmm_energy() oqp.int1e(mol) diff --git a/pyoqp/oqp/library/libdlfind.py b/pyoqp/oqp/library/libdlfind.py index 0800277..45fad44 100644 --- a/pyoqp/oqp/library/libdlfind.py +++ b/pyoqp/oqp/library/libdlfind.py @@ -8,6 +8,7 @@ import numpy as np from oqp.library.libscipy import Optimizer from oqp.utils.file_utils import dump_data, dump_log +import oqp.utils.qmmm as qmmm from libdlfind import dl_find from libdlfind.callback import ( dlf_get_gradient_wrapper, @@ -270,3 +271,98 @@ def check_convergence(self): dump_log(self.mol, title='PyOQP: Geometry Optimization Has Not Converged. Reached The Maximum Iteration') raise StopIteration + +class DLFindQMMM(DLFinder): + """ + OQP DL-Find state specific optimization class + + """ + def __init__(self, mol): + super().__init__(mol) + + dump_log(self.mol, title='PyOQP: DL-FIND Local Minimum QM/MM Optimization', section='dlf') + + self.ref_coord,self.mass,self.natom,self.atoms = qmmm.openmm_info() + self.pre_coord = np.array(self.ref_coord).reshape((-1)) + + self.dlf_get_params = make_dlf_get_params( + coords=self.pre_coord.tolist(), + printl=self.printl, + icoord=self.icoord, + iopt=self.iopt, + imultistate=self.ims, + maxcycle=self.maxit, + tolerance=1e-20, + tolerance_e=1e-20, + ) + + + def one_step(self, coordinates): + # flatten coord + coordinates = coordinates.reshape(-1) + + # add iteration + self.itr += 1 + + dump_log(self.mol, title='PyOQP: QM/MM Geometry Optimization Step %s' % self.itr) + + # update coordinates + qmmm.openmm_update_system(coordinates) + current_xyz=coordinates.reshape((-1,3)) +# coordinates_qm = np.array (()) +# for i in qmmm.qm_atoms: +# coordinates_qm=np.append(coordinates_qm,current_xyz[i][0]) +# coordinates_qm=np.append(coordinates_qm,current_xyz[i][1]) +# coordinates_qm=np.append(coordinates_qm,current_xyz[i][2]) + + num_atoms, x, y, z, q, mass = qmmm.openmm_system() + coordinates_qm=np.array(x + y + z).reshape((3, num_atoms)).T.reshape(-1) + + self.mol.update_system(coordinates_qm) + + if self.itr == 1: + do_init_scf = True + else: + do_init_scf = self.init_scf + oqp.library.ints_1e(self.mol) + + # compute energy + energies = self.sp.energy(do_init_scf=do_init_scf) + + # compute QM/MM gradient + current_xyz = self.mol.get_system().reshape((-1, 3)) + gradient_qm,gradient_mm=qmmm.openmm_gradient(current_xyz,self.mol.data["OQP::partial_charges"]) + self.mol.data["OQP::mm_gradient"]=np.transpose(gradient_qm).tolist() + + self.grad.grads = [self.istate] + grads = self.grad.gradient() + self.mol.energies = energies + self.mol.grads = grads + + qmmm.gradient_qmmm=qmmm.form_gradient_qmmm(grads,gradient_mm) + + # flatten data + energy = energies[self.istate] + grad = qmmm.gradient_qmmm.reshape(-1) + + # evaluate metrics + de = energy - self.pre_energy + rmsd_step = np.mean((coordinates - self.pre_coord) ** 2) ** 0.5 + max_step = np.amax(np.abs(coordinates - self.pre_coord)) + rmsd_grad = np.mean(grad ** 2) ** 0.5 + max_grad = np.amax(np.abs(grad)) + self.metrics['itr'] = self.itr + self.metrics['de'] = de + self.metrics['rmsd_step'] = rmsd_step + self.metrics['max_step'] = max_step + self.metrics['rmsd_grad'] = rmsd_grad + self.metrics['max_grad'] = max_grad + + # store energy and coordinates + self.pre_energy = energy + self.pre_coord = coordinates.copy() + dump_data((self.itr, self.atoms, coordinates, energy, de, rmsd_step, max_step, rmsd_grad, max_grad), + title='OPTIMIZATION', fpath=self.mol.log_path) + + return energy, grad.reshape((self.natom, 3)) + diff --git a/pyoqp/oqp/library/libopenmm.py b/pyoqp/oqp/library/libopenmm.py new file mode 100644 index 0000000..c9ae7ed --- /dev/null +++ b/pyoqp/oqp/library/libopenmm.py @@ -0,0 +1,128 @@ +"""OQP OpenMM MD class""" +import copy +import oqp +import numpy as np +import scipy as sc +from oqp.library.single_point import SinglePoint, Gradient, LastStep +from oqp.utils.file_utils import dump_log, dump_data +import oqp.utils.qmmm as qmmm +import openmm as mm +import openmm.app as app +import openmm.unit as unit + + +class QMMM_MD: + def __init__(self, mol): + self.mol = mol + self.nSteps = mol.config['qmmm']['nsteps'] + self.timeStep = mol.config['qmmm']['timestep'] + self.init_scf = mol.config['optimize']['init_scf'] + self.nstate = mol.config['tdhf']['nstate'] + self.natom = mol.data['natom'] + self.atoms = mol.get_atoms().reshape((self.natom, 1)) + self.sp = SinglePoint(mol) + self.grad = Gradient(mol) + + dump_log(self.mol, title='PyOQP: Entering QM/MM molecular dynamics') + + def run_md(self): + + pdb = qmmm.pdb0 + +#Create an empty system and add particles based on PDB topology + system = mm.System() + for atom in pdb.topology.atoms(): + system.addParticle(atom.element.mass) + number_of_particles = system.getNumParticles() + +#Add QM/MM gradient and energy to system as a CustomExternalForce + oqp_qmmm = mm.CustomExternalForce('grad_x*x + grad_y*y + grad_z*z + qmmm_energy - ecorr') + force_index = system.addForce(oqp_qmmm) + oqp_qmmm.addPerParticleParameter('grad_x') + oqp_qmmm.addPerParticleParameter('grad_y') + oqp_qmmm.addPerParticleParameter('grad_z') + +#First call to OQP and define initial parameters as QM/MM energy and gradient + qmmm_energy, qmmm_grad = self.oqp_qmmm_single_point(0) +# Energy + qmmm_energy = (qmmm_energy/number_of_particles/qmmm.kj_per_mol_to_hartree)*unit.kilojoules_per_mole + oqp_qmmm.addGlobalParameter('qmmm_energy', qmmm_energy) +# Gradient + qmmm_grad = (qmmm_grad/qmmm.kj_per_mol_to_hartree/qmmm.bohr_to_nm) * unit.kilojoules_per_mole/unit.nanometer + ecorr=0.0*unit.kilojoules_per_mole + for i in range(number_of_particles): + oqp_qmmm.addParticle(i, qmmm_grad[i]) + ecorr+=pdb.positions[i][0]*qmmm_grad[i][0] + ecorr+=pdb.positions[i][1]*qmmm_grad[i][1] + ecorr+=pdb.positions[i][2]*qmmm_grad[i][2] + ecorr/=number_of_particles +#This is an energy correction to eliminate the grad_x*x + grad_y*y + grad_z*z contribution + oqp_qmmm.addGlobalParameter('ecorr', ecorr) + +#Create simulation integrator + integrator = mm.VerletIntegrator(self.timeStep*unit.femtoseconds) + simulation = app.Simulation(pdb.topology, system, integrator) + simulation.context.setPositions(pdb.positions) + simulation.context.setVelocitiesToTemperature(300.0*unit.kelvin) + +#Add reporters + simulation.reporters.append(app.PDBReporter('trajectory.pdb', 1)) + oqp_output=open(self.mol.log+'.md','w') + oqp_output.write(f"\n\nStarting NVE dynamics using Verlet algorithm:\n") + simulation.reporters.append(app.StateDataReporter(oqp_output, 1, step=True, time=True, separator=" ",totalEnergy=True, kineticEnergy=True, potentialEnergy=True, temperature=False, speed=True,elapsedTime=True)) + + state = simulation.context.getState(getPositions=True,getForces=True) + + for step in range(self.nSteps): + + simulation.step(1) + +#Get the current positions from simulation context and update the pdb + state = simulation.context.getState(getPositions=True) + qmmm.pdb0.positions = state.getPositions()[:number_of_particles] + +#Perform OQP calculation and update parameters in context + qmmm_energy, qmmm_grad = self.oqp_qmmm_single_point(step+1) + qmmm_energy = qmmm_energy/number_of_particles*unit.kilojoules_per_mole/qmmm.kj_per_mol_to_hartree + simulation.context.setParameter('qmmm_energy', qmmm_energy) + + qmmm_grad = (qmmm_grad/qmmm.kj_per_mol_to_hartree/qmmm.bohr_to_nm) * unit.kilojoules_per_mole/unit.nanometer + ecorr=0.0*unit.kilojoules_per_mole + for i in range(number_of_particles): + oqp_qmmm.setParticleParameters(i, i, qmmm_grad[i]) + ecorr+=pdb.positions[i][0]*qmmm_grad[i][0] + ecorr+=pdb.positions[i][1]*qmmm_grad[i][1] + ecorr+=pdb.positions[i][2]*qmmm_grad[i][2] + ecorr/=number_of_particles + simulation.context.setParameter('ecorr', ecorr) + oqp_qmmm.updateParametersInContext(simulation.context) + + def oqp_qmmm_single_point(self,itr): + + dump_log(self.mol, title='PyOQP: QM/MM molecular dynamics (NVE) using Verlet algorithm ') + + num_atoms, x, y, z, q, mass = qmmm.openmm_system() + coordinates_qm=np.array(x + y + z).reshape((3, num_atoms)).T.reshape(-1) + + self.mol.update_system(coordinates_qm) + + # compute 1e integral + if itr == 0: + do_init_scf = True + else: + do_init_scf = self.init_scf + oqp.library.ints_1e(self.mol) + + # compute energy + energies = self.sp.energy(do_init_scf=do_init_scf) + + # compute QM/MM gradient + grads = self.grad.gradient() + + # flatten data + energy = energies[qmmm.istate] + grad = qmmm.gradient_qmmm + + return energy, grad + + diff --git a/pyoqp/oqp/library/libscipy.py b/pyoqp/oqp/library/libscipy.py index 32b782f..80c1a38 100644 --- a/pyoqp/oqp/library/libscipy.py +++ b/pyoqp/oqp/library/libscipy.py @@ -5,7 +5,7 @@ import scipy as sc from oqp.library.single_point import SinglePoint, Gradient, LastStep from oqp.utils.file_utils import dump_log, dump_data - +import oqp.utils.qmmm as qmmm class Optimizer: """ @@ -925,3 +925,124 @@ def ReLU(x, y): return x else: return 0 + +class QMMMOpt(Optimizer): + """ + OQP QM/MM optimization class + + """ + + def __init__(self, mol): + super().__init__(mol) + self.mol = mol + self.metrics['radius'] = 0 + self.metrics['distance'] = 0 + self.metrics['step_size'] = self.step_size + self.metrics['step_tol'] = self.step_tol + self.ref_coord,self.mass,self.natom,self.atoms = qmmm.openmm_info() + self.pre_coord = self.ref_coord + self.tmass = np.sum(self.mass) + self.init_energy = 0 + self.last_energy = 0 + self.init_xyz = None + self.last_xyz = None + self.last_radius = 0 + self.opt_status = 0 + self.message = 'not converged' + + if not self.mol.config['input']['qmmm_flag']: + exit(f"QM/MM optimizations require an active qmmm_flag") + + def one_step(self, coordinates): + # add iteration + self.itr += 1 + + dump_log(self.mol, title='PyOQP: QM/MM Geometry Optimization Step %s' % self.itr) + + + # update coordinates + qmmm.openmm_update_system(coordinates) + current_xyz=coordinates.reshape((-1,3)) +# for i in qmmm.qm_atoms: +# coordinates_qm=np.append(coordinates_qm,current_xyz[i][0]) +# coordinates_qm=np.append(coordinates_qm,current_xyz[i][1]) +# coordinates_qm=np.append(coordinates_qm,current_xyz[i][2]) + + num_atoms, x, y, z, q, mass = qmmm.openmm_system() + coordinates_qm=np.array(x + y + z).reshape((3, num_atoms)).T.reshape(-1) + + self.mol.update_system(coordinates_qm) + + # compute 1e integral + if self.itr == 1: + do_init_scf = True + else: + do_init_scf = self.init_scf + oqp.library.ints_1e(self.mol) + + # compute energy + energies = self.sp.energy(do_init_scf=do_init_scf) + + # compute QM/MM gradient + current_xyz = self.mol.get_system().reshape((-1, 3)) + gradient_qm,gradient_mm=qmmm.openmm_gradient(current_xyz,self.mol.data["OQP::partial_charges"]) + self.mol.data["OQP::mm_gradient"]=np.transpose(gradient_qm).tolist() + + self.grad.grads = [self.istate] + grads = self.grad.gradient() + self.mol.energies = energies + self.mol.grads = grads + + qmmm.gradient_qmmm=qmmm.form_gradient_qmmm(grads,gradient_mm) + + # flatten data + energy = energies[self.istate] + grad = qmmm.gradient_qmmm.reshape(-1) + + # evaluate metrics + de = energy - self.pre_energy + rmsd_step = np.mean((coordinates - self.pre_coord) ** 2) ** 0.5 + max_step = np.amax(np.abs(coordinates - self.pre_coord)) + rmsd_grad = np.mean(grad ** 2) ** 0.5 + max_grad = np.amax(np.abs(grad)) + self.metrics['itr'] = self.itr + self.metrics['de'] = de + self.metrics['rmsd_step'] = rmsd_step + self.metrics['max_step'] = max_step + self.metrics['rmsd_grad'] = rmsd_grad + self.metrics['max_grad'] = max_grad + + # store energy and coordinates + self.pre_energy = energy + self.pre_coord = coordinates.copy() + dump_data((self.itr, self.atoms, coordinates, energy, de, rmsd_step, max_step, rmsd_grad, max_grad), + title='QM/MM OPTIMIZATION', fpath=self.mol.log_path) + + return energy, grad + + def check_convergence(self): + # write convergence to log + dump_log( + self.mol, title='Geometry Optimization Convergence %s' % self.itr, + section='QM/MM', info=self.metrics + ) + + if np.abs(self.metrics['de']) <= self.energy_shift and \ + self.metrics['rmsd_step'] <= self.rmsd_step and \ + self.metrics['max_step'] <= self.rmsd_step and \ + self.metrics['rmsd_grad'] <= self.rmsd_grad and \ + self.metrics['max_grad'] <= self.max_grad and \ + self.metrics['radius'] > self.step_tol: + dump_log(self.mol, title='PyOQP: Geometry Optimization Has Converged') + self.opt_status = 1 + self.message = 'converged' + raise StopIteration + + else: + if self.itr == self.maxit: + dump_log(self.mol, + title='PyOQP: Geometry Optimization Has Not Converged. Reached The Maximum Iteration') + raise StopIteration + + + diff --git a/pyoqp/oqp/library/runfunc.py b/pyoqp/oqp/library/runfunc.py index cae7909..d6c65ea 100644 --- a/pyoqp/oqp/library/runfunc.py +++ b/pyoqp/oqp/library/runfunc.py @@ -7,8 +7,9 @@ BasisOverlap, NACME, NAC ) -from oqp.library.libscipy import StateSpecificOpt, MECIOpt, MECPOpt, MEP -from oqp.library.libdlfind import DLFindMin, DLFindTS, DLFindMECI +from oqp.library.libscipy import StateSpecificOpt, MECIOpt, MEP, QMMMOpt +from oqp.library.libdlfind import DLFindMin, DLFindTS, DLFindMECI, DLFindQMMM +from oqp.library.libopenmm import QMMM_MD def prep_guess(mol): @@ -45,6 +46,10 @@ def compute_scf_prop(mol): else: raise ValueError(f'Unknown property: {prop}') + if 'resp' not in properties and mol.config['input']['qmmm_flag']: + oqp.resp_charges(mol) + + def compute_grad(mol): # prepare guess orbital @@ -56,9 +61,22 @@ def compute_grad(mol): # compute gradient Gradient(mol).gradient() + # compute properties + if mol.config['input']['qmmm_flag']: + oqp.resp_charges(mol) + # compute dftd4 LastStep(mol).compute(mol, grad_list=mol.config['properties']['grad']) +def compute_md(mol): + + # prepare guess orbital + prep_guess(mol) + + #Run MD + qmmm_md = QMMM_MD(mol) + qmmm_md.run_md() + def compute_nacme(mol): # prepare guess orbital prep_guess(mol) @@ -213,6 +231,10 @@ def get_optimizer(mol): }, } + if runtype == 'optimize' and mol.config['input']['qmmm_flag']: + if lib == 'dlfind': opt_lib[lib][runtype] = DLFindQMMM + elif lib == 'scipy': opt_lib[lib][runtype] = QMMMOpt + if opt_lib[lib][runtype]: return opt_lib[lib][runtype](mol) diff --git a/pyoqp/oqp/library/single_point.py b/pyoqp/oqp/library/single_point.py index e730379..9589f69 100644 --- a/pyoqp/oqp/library/single_point.py +++ b/pyoqp/oqp/library/single_point.py @@ -21,7 +21,7 @@ from oqp.library.frequency import normal_mode, thermal_analysis from oqp.utils.file_utils import dump_log, dump_data, write_config, write_xyz - +import oqp.utils.qmmm as qmmm class Calculator: """ @@ -390,6 +390,11 @@ def gradient(self): dump_log(self.mol, title='PyOQP: Entering Gradient Calculation') + if self.mol.config['input']['qmmm_flag']: + current_xyz = self.mol.get_system().reshape((-1, 3)) + gradient_qm,gradient_mm=qmmm.openmm_gradient(current_xyz,self.mol.data["OQP::partial_charges"]) + self.mol.data["OQP::mm_gradient"]=gradient_qm + # compute gradients grads = [] if self.method == 'hf': @@ -399,6 +404,9 @@ def gradient(self): self.mol.grads = grads + if self.mol.config['input']['qmmm_flag']: + qmmm.gradient_qmmm=qmmm.form_gradient_qmmm(grads,gradient_mm) + return grads def scf_grad(self): diff --git a/pyoqp/oqp/molecule/molecule.py b/pyoqp/oqp/molecule/molecule.py index 5f6112e..876646d 100644 --- a/pyoqp/oqp/molecule/molecule.py +++ b/pyoqp/oqp/molecule/molecule.py @@ -59,6 +59,7 @@ def __init__(self, project_name, input_file, log, 'OQP::td_abxc', 'OQP::td_bvec_mo', 'OQP::td_mrsf_density', 'OQP::td_energies', 'OQP::td_states_overlap', 'OQP::dc_matrix', 'OQP::nac_matrix', + 'OQP::hamiltonian_qmmm', 'OQP::mm_potential', 'OQP::charge_operator', 'OQP::partial_charges' ] self.start_time = None diff --git a/pyoqp/oqp/molecule/oqpdata.py b/pyoqp/oqp/molecule/oqpdata.py index 7064dfa..7ec73d9 100644 --- a/pyoqp/oqp/molecule/oqpdata.py +++ b/pyoqp/oqp/molecule/oqpdata.py @@ -5,6 +5,7 @@ from oqp import ffi, lib from oqp.periodic_table import MASSES, SYMBOL_MAP from oqp.utils.constants import ANGSTROM_TO_BOHR +import oqp.utils.qmmm as qmmm def sarray(strng): @@ -44,6 +45,15 @@ def path(strng): OQP_CONFIG_SCHEMA = { + 'qmmm': { + 'forcefield': {'type': sarray, 'default': 'amber14-all.xml,amber14/tip3p.xml'}, + 'nonbondedmethod': {'type': str, 'default': 'NoCutoff'}, + 'constraints': {'type': str, 'default': 'None'}, + 'rigidwater': {'type': bool, 'default': 'False'}, + 'nsteps': {'type': int, 'default': '1'}, + 'timestep': {'type': int, 'default': '1'}, + 'istate': {'type': int, 'default': '0'}, + }, 'input': { 'charge': {'type': int, 'default': '0'}, 'basis': {'type': string, 'default': ''}, @@ -53,6 +63,7 @@ def path(strng): 'system': {'type': str, 'default': ''}, 'system2': {'type': str, 'default': ''}, 'd4': {'type': bool, 'default': 'False'}, + 'qmmm_flag': {'type': bool, 'default': 'False'}, }, 'guess': { 'type': {'type': string, 'default': 'huckel'}, @@ -205,6 +216,7 @@ class OQPData: "functional": "set_dft_functional", "system": "set_system", "system2": "set_system2", + "qmmm_flag": "set_qmmm_flag", }, "guess": { }, @@ -260,6 +272,15 @@ class OQPData: "spc_coov": "set_tdhf_spc_coov", "conf_threshold": "set_conf_threshold", }, + "qmmm": { + "forcefield": "set_qmmm_forcefield", + "nonbondedmethod": "set_qmmm_nonbondedmethod", + "constraints": "set_qmmm_constraints", + "rigidwater": "set_qmmm_rigidwater", + "nsteps": "set_qmmm_nsteps", + "timestep": "set_qmmm_timestep", + "istate": "set_qmmm_istate", + }, } _typemap = [np.void, np.int8, @@ -449,6 +470,38 @@ def set_scf_incremental(self, flag): """Set incremental Fock matrix build""" self._data.control.scf_incremental = 1 if flag else 0 + def set_qmmm_flag(self, qmmm_flag): + """Handle QM/MM calculation type""" + self._data.control.qmmm_flag=qmmm_flag + + def set_qmmm_forcefield(self, forcefield): + """Handle QM/MM calculation forcefield""" + qmmm.force_field = forcefield + + def set_qmmm_rigidwater(self, rigidwater): + """Handle QM/MM calculation rigidWater""" + qmmm.rigidWater = rigidwater + + def set_qmmm_nonbondedmethod(self, nonbondedmethod): + """Handle QM/MM calculation nonbondedMethod""" + qmmm.nonbondedMethod = nonbondedmethod + + def set_qmmm_constraints(self, constraints): + """Handle QM/MM calculation constraints""" + qmmm.contraints = constraints + + def set_qmmm_nsteps(self, nsteps): + """Handle QM/MM calculation constraints""" + qmmm.nSteps = nsteps + + def set_qmmm_timestep(self, timestep): + """Handle QM/MM calculation constraints""" + qmmm.timeStep = timestep + + def set_qmmm_istate(self, istate): + """Handle QM/MM calculation istate""" + qmmm.istate = istate + def set_tdhf_type(self, td_type): """Handle td-dft calculation type""" if td_type.lower() == 'tda': @@ -699,8 +752,10 @@ def compute_alpha_beta_electrons(n_e, mult): def read_system(system): - system = system.split("\n") - if system[0]: + system0 = system + system = system.split() + """Set up atomic data""" + if system[0].lower().endswith('.xyz'): if not os.path.exists(system[0]): raise FileNotFoundError("XYZ file %s is not found!" % system[0]) @@ -709,22 +764,57 @@ def read_system(system): num_atoms = int(system[0]) system = system[2: 2 + num_atoms] + atoms = [] + for i, line in enumerate(system): + line = line.split() + if len(line) >= 4: + atoms.append(line[0: 4]) + else: + print(f"{system[i]} is not valid line for atom configuration!") + + q = [float(SYMBOL_MAP[atoms[i][0]]) for i in range(0, num_atoms)] + x = [float(atoms[i][1]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] + y = [float(atoms[i][2]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] + z = [float(atoms[i][3]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] + mass = [MASSES[int(SYMBOL_MAP[atoms[i][0]])] for i in range(0, num_atoms)] + elif system[0].lower().endswith('.pdb'): + + if not os.path.exists(system[0]): + raise FileNotFoundError("PDB file %s is not found!" % system[0]) + qmmm.pdb_file=system[0] + + atom_list = [] + for i in system[1:]: + if i.find('-') != -1: + start, end = map(int, i.split('-')) + atom_list.extend(list(range(start, end + 1))) + else: + atom_list.append(int(i)) + + if len(atom_list) != len(set(atom_list)): + raise ValueError("Repeated entries in QM atom list") + + if any(value < 0 for value in atom_list): + raise ValueError("Negative indexes are not allowed in QM atom list") + + qmmm.qm_atoms,qmmm.pdb0,qmmm.forcefield0,qmmm.system0=qmmm.openmm_init(atom_list=atom_list) + num_atoms, x, y, z, q, mass = qmmm.openmm_system() else: + system = system0.split("\n") system = system[1:] num_atoms = len(system) + atoms = [] + for i, line in enumerate(system): + line = line.split() + if len(line) >= 4: + atoms.append(line[0: 4]) + else: + print(f"{system[i]} is not valid line for atom configuration!") - atoms = [] - for i, line in enumerate(system): - line = line.split() - if len(line) >= 4: - atoms.append(line[0: 4]) - else: - print(f"{system[i]} is not valid line for atom configuration!") - - q = [float(SYMBOL_MAP[atoms[i][0]]) for i in range(0, num_atoms)] - x = [float(atoms[i][1]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] - y = [float(atoms[i][2]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] - z = [float(atoms[i][3]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] - mass = [MASSES[int(SYMBOL_MAP[atoms[i][0]])] for i in range(0, num_atoms)] + q = [float(SYMBOL_MAP[atoms[i][0]]) for i in range(0, num_atoms)] + x = [float(atoms[i][1]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] + y = [float(atoms[i][2]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] + z = [float(atoms[i][3]) / ANGSTROM_TO_BOHR for i in range(0, num_atoms)] + mass = [MASSES[int(SYMBOL_MAP[atoms[i][0]])] for i in range(0, num_atoms)] return num_atoms, x, y, z, q, mass diff --git a/pyoqp/oqp/pyoqp.py b/pyoqp/oqp/pyoqp.py index 81cd5ab..32a2d56 100644 --- a/pyoqp/oqp/pyoqp.py +++ b/pyoqp/oqp/pyoqp.py @@ -14,7 +14,7 @@ from oqp.utils.input_checker import check_input_values from oqp.molecule import Molecule from oqp.library.runfunc import ( - compute_energy, compute_grad, compute_nac, compute_soc, compute_geom, + compute_energy, compute_grad, compute_nac, compute_soc, compute_geom, compute_md, compute_nacme, compute_properties, compute_data, compute_hess, compute_thermo ) from oqp.utils.mpi_utils import MPIManager @@ -61,6 +61,7 @@ def __init__(self, project=None, input_file=None, log=None, 'meci': compute_geom, 'mecp': compute_geom, 'mep': compute_geom, + 'md': compute_md, 'ts': compute_geom, 'hess': compute_hess, 'thermo': compute_thermo, diff --git a/pyoqp/oqp/utils/file_utils.py b/pyoqp/oqp/utils/file_utils.py index 6d45fa4..dd896c4 100644 --- a/pyoqp/oqp/utils/file_utils.py +++ b/pyoqp/oqp/utils/file_utils.py @@ -8,7 +8,7 @@ from oqp.periodic_table import SYMBOL_MAP, ELEMENTS_NAME from oqp.utils.constants import ANGSTROM_TO_BOHR from oqp.utils.mpi_utils import mpi_dump - +from oqp.utils.qmmm import gradient_qmmm def try_basis(basis, path=None, fallback='6-31g'): """try various basis file locations and return the matching one""" @@ -225,19 +225,20 @@ def dump_log(mol, title=None, section=None, info=None, must_print=False): loginfo += f' PyOQP state {n:<6} {energy:<16.8f}\n' if section == 'grad': - atoms = mol.get_atoms() - loginfo += ' PyOQP electronic gradients\n' - for n in info['grad_list']: - grad = write_grad(atoms, info['el'][n]) - loginfo += f' PyOQP state {n:<6}\n{grad}\n' - - d4 = info['d4'] - loginfo += f'\n PyOQP dftd correction\n' - loginfo += f'{write_grad(atoms, d4)}\n\n' - loginfo += ' PyOQP dispersion corrected gradients\n' - for n in info['grad_list']: - grad = write_grad(atoms, mol.grads[n]) - loginfo += f' PyOQP state {n:<6}\n{grad}\n' + if not mol.config['input']['qmmm_flag']: + atoms = mol.get_atoms() + loginfo += ' PyOQP electronic gradients\n' + for n in info['grad_list']: + grad = write_grad(atoms, info['el'][n]) + loginfo += f' PyOQP state {n:<6}\n{grad}\n' + + d4 = info['d4'] + loginfo += f'\n PyOQP dftd correction\n' + loginfo += f'{write_grad(atoms, d4)}\n\n' + loginfo += ' PyOQP dispersion corrected gradients\n' + for n in info['grad_list']: + grad = write_grad(atoms, mol.grads[n]) + loginfo += f' PyOQP state {n:<6}\n{grad}\n' if section == 'opt': loginfo += """ diff --git a/pyoqp/oqp/utils/input_checker.py b/pyoqp/oqp/utils/input_checker.py index 6de04d0..58f49fa 100644 --- a/pyoqp/oqp/utils/input_checker.py +++ b/pyoqp/oqp/utils/input_checker.py @@ -1,265 +1,276 @@ -"""Input value checker""" - -import os -import multiprocessing -from oqp.utils.mpi_utils import MPIManager - - -def check_input_values(config): - runtype = config['input']['runtype'] - - check_func = { - 'energy': check_energy_input, - 'grad': check_grad_input, - 'nac': check_nac_input, - 'bp': check_nac_input, - 'soc': check_soc_input, - 'optimize': check_optimize_input, - 'meci': check_optimize_input, - 'mecp': check_optimize_input, - 'mep': check_optimize_input, - 'ts': check_optimize_input, - 'neb': check_optimize_input, - 'hess': check_hess_input, - 'nacme': check_nacme_input, - 'prop': skip_check, - 'data': skip_check, - } - - info = f'\nPyOQP checking input\n\n[input] runtype={runtype}' - check_func[runtype](config, info) - - -def skip_check(config, info): - pass - -def not_available(func_name, info): - exit(f'{info}\nPyOQP: runtype {func_name} is not available yet\n') - -def check_scf_input(config, info): - scf_type = config['scf']['type'] - scf_mult = config['scf']['multiplicity'] - - info += f'\n[scf] type={scf_type}' - info += f'\n[scf] multiplicity={scf_mult}' - - if scf_mult < 1: - exit(f'{info}\nPyOQP:scf multiplicity must be greater than 0\n') - - if scf_mult > 1 and scf_type == 'rhf': - exit(f'{info}\nPyOQP: scf multiplicity {scf_mult} does not match rhf, choose uhf or rohf\n') - -def check_tdhf_input(config, info): - check_scf_input(config, info) - - scf_type = config['scf']['type'] - scf_mult = config['scf']['multiplicity'] - td_type = config['tdhf']['type'] - td_mult = config['tdhf']['multiplicity'] - - info += f'\n[scf] type={scf_type}' - info += f'\n[scf] multiplicity={scf_mult}' - info += f'\n[tdhf] type={td_type}' - info += f'\n[tdhf] multiplicity={td_mult}' - - if td_type in ['rpa', 'tda'] and scf_mult != td_mult: - print(f'{info}\nPyOQP: Caution! tdhf type {td_type} multiplicity {td_mult} is not equal to scf multiplicity {scf_mult}\n') - - if td_type in ['sf', 'mrsf'] and scf_mult == td_mult: - print(f'{info}\nPyOQP: Caution! tdhf type {td_type} multiplicity {td_mult} is equal to scf multiplicity {scf_mult}\n') - - if td_type in ['mrsf', 'sf'] and scf_type != 'rohf': - exit(f'{info}\nPyOQP: tdhf type {td_type} cannot use {scf_type} orbitals, choose rohf in scf type\n') - -def check_energy_input(config, info): - method = config['input']['method'] - - check_func = { - 'hf': check_scf_input, - 'tdhf': check_tdhf_input, - } - - info += f'\n[input] method={method}' - check_func[method](config, info) - -def check_grad_input(config, info): - - check_energy_input(config, info) - - method = config['input']['method'] - td_type = config['tdhf']['type'] - grad = config['properties']['grad'] - - info += f'\n[input] method={method}' - info += f'\n[tdhf] type={td_type}' - info += f'\n[properties] grad={grad}' - - if method == 'hf' and sum(grad) > 0: - exit(f'{info}\nPyOQP: scf cannot compute gradient for state > 0\n') - - if method == 'tdhf' and 0 in grad: - exit(f'{info}\nPyOQP: tdhf type {td_type} cannot compute gradient for state 0 \n') - -def check_optimize_input(config, info): - runtype = config['input']['runtype'] - check_energy_input(config, info) - method = config['input']['method'] - istate = config['optimize']['istate'] - jstate = config['optimize']['jstate'] - imult = config['optimize']['imult'] - jmult = config['optimize']['jmult'] - - info += f'\n[input] method={method}' - info += f'\n[optimize] istate={istate}' - - if method == 'hf' and istate > 0: - exit(f'{info}\nPyOQP: scf cannot compute gradient for istate > 0\n') - - if method == 'tdhf' and istate == 0: - exit(f'{info}\nPyOQP: tdhf cannot compute gradient for istate = 0 \n') - - info += f'\n[optimize] jstate={jstate}' - if runtype == 'meci' and jstate - istate < 1: - exit(f'{info}\nPyOQP: meci optimization require jstate >= istate + 1 \n') - - if runtype == 'mecp' and imult == jmult: - exit(f'{info}\nPyOQP: mecp optimization require imult != jmult \n') - - lib = config['optimize']['lib'] - info += f'\n[optimize] lib={lib}' - - if lib not in ['scipy', 'dlfind']: - exit(f'{info}\nPyOQP: cannot recognize the geometry optimization library {lib} \n') - - if lib == 'dlfind': - check_dlfind_input(config, info) - -def check_dlfind_input(config, info): - runtype = config['input']['runtype'] - icoord = config['dlfind']['icoord'] - iopt = config['dlfind']['iopt'] - ims = config['dlfind']['ims'] - - info += f'\n[dlfind] icoord={icoord}' - info += f'\n[dlfind] iopt={iopt}' - info += f'\n[dlfind] ims={ims}' - - check_func = { - 'optimize': check_dlfind_min_input, - 'meci': check_dlfind_meci_input, - 'ts': check_dlfind_ts_input, - } - - check_func[runtype](icoord, iopt, ims, info) - -def check_dlfind_min_input(icoord, iopt, ims, info): - if icoord not in [0, 1, 2, 3, 4]: - exit(f'{info}\nPyOQP: dl-find single state optimization only support icoord 0–4, but found {icoord}') - - if iopt not in [0, 1, 2, 3]: - exit(f'{info}\nPyOQP: dl-find single state optimization only support iopt 0–3, but found {iopt}') - - if ims not in [0]: - exit(f'{info}\nPyOQP: dl-find single state optimization only support ims 0, but found {ims}') - -def check_dlfind_meci_input(icoord, iopt, ims, info): - if iopt not in [0, 1, 2, 3]: - exit(f'{info}\nPyOQP: dl-find meci optimization only support iopt 0–3, but found {iopt}') - - if ims not in [1, 2, 3]: - exit(f'{info}\nPyOQP: dl-find meci optimization only support ims 1-3, but found {ims}') - - if ims == 3 and icoord not in [10, 11, 12, 13, 14]: - exit(f'{info}\nPyOQP: dl-find meci optimization {ims} only support icoord 10–14, but found {icoord}') - - if ims in [1, 2] and icoord not in [0, 1, 2, 3, 4]: - exit(f'{info}\nPyOQP: dl-find meci optimization {ims} only support icoord 0–4, but found {icoord}') - -def check_dlfind_ts_input(icoord, iopt, ims, info): - if iopt not in [0, 1, 2, 3]: - exit(f'{info}\nPyOQP: dl-find meci optimization only support iopt 0–3, but found {iopt}') - - if ims not in [1, 2, 3]: - exit(f'{info}\nPyOQP: dl-find meci optimization only support ims 1-3, but found {ims}') - - if ims == 3 and icoord not in [10, 11, 12, 13, 14]: - exit(f'{info}\nPyOQP: dl-find meci optimization {ims} only support icoord 10–14, but found {icoord}') - - if ims in [1, 2] and icoord not in [0, 1, 2, 3, 4]: - exit(f'{info}\nPyOQP: dl-find meci optimization {ims} only support icoord 0–4, but found {icoord}') - - -def check_nac_input(config, info): - check_energy_input(config, info) - check_energy_input(config, info) - method = config['input']['method'] - td = config['tdhf']['type'] - nproc = config['nac']['nproc'] - ncpu = multiprocessing.cpu_count() - mpi = MPIManager().size - omp = int(os.environ['OMP_NUM_THREADS']) - - info += f'\n[input] method={method}' - - if method == 'hf': - exit(f'{info}\nPyOQP: scf cannot compute nac for hf calculations\n') - - if method == 'tdhf' and td != 'mrsf': - exit(f'{info}\nPyOQP: tdhf cannot compute nac for non-mrsf calculations\n') - - if MPIManager().use_mpi: - info += f'\nMPI process={mpi}' - info += f'\nOMP_NUM_THREADS={omp}' - print(f'{info}\nPyOQP: nac requested {mpi * omp} cpus! {ncpu} available on a single node') - else: - info += f'\n[nac] nproc={nproc}' - info += f'\nOMP_NUM_THREADS={omp}' - print(f'{info}\nPyOQP: nac requested {nproc * omp} cpus! {ncpu} available') - - if nproc < 1: - exit(f'{info}\nPyOQP: nac requested {nproc} process is illegal') - - -def check_soc_input(config, info): - not_available('soc', info) - - -def check_neb_input(config, info): - not_available('neb', info) - - -def check_hess_input(config, info): - check_energy_input(config, info) - method = config['input']['method'] - state = config['hess']['state'] - nproc = config['hess']['nproc'] - restart = config['hess']['restart'] - ncpu = multiprocessing.cpu_count() - mpi = MPIManager().size - omp = int(os.environ['OMP_NUM_THREADS']) - - info += f'\n[input] method={method}' - info += f'\n[hess] state={state}' - - if method == 'hf' and state > 0: - exit(f'{info}\nPyOQP: scf cannot compute hessian for state > 0\n') - - if method == 'tdhf' and state == 0: - exit(f'{info}\nPyOQP: tdhf cannot compute hessian for state = 0 \n') - - if restart != 'read': - if MPIManager().use_mpi: - info += f'\nMPI process={mpi}' - info += f'\nOMP_NUM_THREADS={omp}' - print(f'{info}\nPyOQP: hessian requested {mpi * omp} cpus! {ncpu} available on a single node') - else: - info += f'\n[hess] nproc={nproc}' - info += f'\nOMP_NUM_THREADS={omp}' - print(f'{info}\nPyOQP: hessian requested {nproc * omp} cpus! {ncpu} available') - - if nproc < 1: - exit(f'{info}\nPyOQP: hessian requested {nproc} process is illegal') - - -def check_nacme_input(config, info): - check_energy_input(config, info) +"""Input value checker""" + +import os +import multiprocessing +import oqp.utils.qmmm as qmmm +from oqp.utils.mpi_utils import MPIManager + + +def check_input_values(config): + runtype = config['input']['runtype'] + + check_func = { + 'energy': check_energy_input, + 'grad': check_grad_input, + 'nac': check_nac_input, + 'bp': check_nac_input, + 'soc': check_soc_input, + 'optimize': check_optimize_input, + 'md': check_md_input, + 'meci': check_optimize_input, + 'mecp': check_optimize_input, + 'mep': check_optimize_input, + 'ts': check_optimize_input, + 'neb': check_optimize_input, + 'hess': check_hess_input, + 'nacme': check_nacme_input, + 'qmmm': check_qmmm_input, + 'prop': skip_check, + 'data': skip_check, + } + + info = f'\nPyOQP checking input\n\n[input] runtype={runtype}' + check_func[runtype](config, info) + + +def skip_check(config, info): + pass + +def not_available(func_name, info): + exit(f'{info}\nPyOQP: runtype {func_name} is not available yet\n') + +def check_md_input(config, info): + if not config['input']['qmmm_flag']: exit(f"MD is only available for QM/MM") + +def check_scf_input(config, info): + scf_type = config['scf']['type'] + scf_mult = config['scf']['multiplicity'] + + info += f'\n[scf] type={scf_type}' + info += f'\n[scf] multiplicity={scf_mult}' + + if scf_mult < 1: + exit(f'{info}\nPyOQP:scf multiplicity must be greater than 0\n') + + if scf_mult > 1 and scf_type == 'rhf': + exit(f'{info}\nPyOQP: scf multiplicity {scf_mult} does not match rhf, choose uhf or rohf\n') + +def check_tdhf_input(config, info): + check_scf_input(config, info) + + scf_type = config['scf']['type'] + scf_mult = config['scf']['multiplicity'] + td_type = config['tdhf']['type'] + td_mult = config['tdhf']['multiplicity'] + + info += f'\n[scf] type={scf_type}' + info += f'\n[scf] multiplicity={scf_mult}' + info += f'\n[tdhf] type={td_type}' + info += f'\n[tdhf] multiplicity={td_mult}' + + if td_type in ['rpa', 'tda'] and scf_mult != td_mult: + print(f'{info}\nPyOQP: Caution! tdhf type {td_type} multiplicity {td_mult} is not equal to scf multiplicity {scf_mult}\n') + + if td_type in ['sf', 'mrsf'] and scf_mult == td_mult: + print(f'{info}\nPyOQP: Caution! tdhf type {td_type} multiplicity {td_mult} is equal to scf multiplicity {scf_mult}\n') + + if td_type in ['mrsf', 'sf'] and scf_type != 'rohf': + exit(f'{info}\nPyOQP: tdhf type {td_type} cannot use {scf_type} orbitals, choose rohf in scf type\n') + +def check_energy_input(config, info): + method = config['input']['method'] + + check_func = { + 'hf': check_scf_input, + 'tdhf': check_tdhf_input, + } + + info += f'\n[input] method={method}' + check_func[method](config, info) + +def check_grad_input(config, info): + + check_energy_input(config, info) + + method = config['input']['method'] + td_type = config['tdhf']['type'] + grad = config['properties']['grad'] + + info += f'\n[input] method={method}' + info += f'\n[tdhf] type={td_type}' + info += f'\n[properties] grad={grad}' + + if method == 'hf' and sum(grad) > 0: + exit(f'{info}\nPyOQP: scf cannot compute gradient for state > 0\n') + + if method == 'tdhf' and 0 in grad: + exit(f'{info}\nPyOQP: tdhf type {td_type} cannot compute gradient for state 0 \n') + +def check_optimize_input(config, info): + runtype = config['input']['runtype'] + check_energy_input(config, info) + method = config['input']['method'] + istate = config['optimize']['istate'] + jstate = config['optimize']['jstate'] + imult = config['optimize']['imult'] + jmult = config['optimize']['jmult'] + + info += f'\n[input] method={method}' + info += f'\n[optimize] istate={istate}' + + if method == 'hf' and istate > 0: + exit(f'{info}\nPyOQP: scf cannot compute gradient for istate > 0\n') + + if method == 'tdhf' and istate == 0: + exit(f'{info}\nPyOQP: tdhf cannot compute gradient for istate = 0 \n') + + info += f'\n[optimize] jstate={jstate}' + if runtype == 'meci' and jstate - istate < 1: + exit(f'{info}\nPyOQP: meci optimization require jstate >= istate + 1 \n') + + if runtype == 'mecp' and imult == jmult: + exit(f'{info}\nPyOQP: mecp optimization require imult != jmult \n') + + lib = config['optimize']['lib'] + info += f'\n[optimize] lib={lib}' + + if lib not in ['scipy', 'dlfind']: + exit(f'{info}\nPyOQP: cannot recognize the geometry optimization library {lib} \n') + + if lib == 'dlfind': + check_dlfind_input(config, info) + +def check_dlfind_input(config, info): + runtype = config['input']['runtype'] + icoord = config['dlfind']['icoord'] + iopt = config['dlfind']['iopt'] + ims = config['dlfind']['ims'] + + info += f'\n[dlfind] icoord={icoord}' + info += f'\n[dlfind] iopt={iopt}' + info += f'\n[dlfind] ims={ims}' + + check_func = { + 'optimize': check_dlfind_min_input, + 'meci': check_dlfind_meci_input, + 'ts': check_dlfind_ts_input, + } + + check_func[runtype](icoord, iopt, ims, info) + +def check_dlfind_min_input(icoord, iopt, ims, info): + if icoord not in [0, 1, 2, 3, 4]: + exit(f'{info}\nPyOQP: dl-find single state optimization only support icoord 0–4, but found {icoord}') + + if iopt not in [0, 1, 2, 3]: + exit(f'{info}\nPyOQP: dl-find single state optimization only support iopt 0–3, but found {iopt}') + + if ims not in [0]: + exit(f'{info}\nPyOQP: dl-find single state optimization only support ims 0, but found {ims}') + +def check_dlfind_meci_input(icoord, iopt, ims, info): + if iopt not in [0, 1, 2, 3]: + exit(f'{info}\nPyOQP: dl-find meci optimization only support iopt 0–3, but found {iopt}') + + if ims not in [1, 2, 3]: + exit(f'{info}\nPyOQP: dl-find meci optimization only support ims 1-3, but found {ims}') + + if ims == 3 and icoord not in [10, 11, 12, 13, 14]: + exit(f'{info}\nPyOQP: dl-find meci optimization {ims} only support icoord 10–14, but found {icoord}') + + if ims in [1, 2] and icoord not in [0, 1, 2, 3, 4]: + exit(f'{info}\nPyOQP: dl-find meci optimization {ims} only support icoord 0–4, but found {icoord}') + +def check_dlfind_ts_input(icoord, iopt, ims, info): + if iopt not in [0, 1, 2, 3]: + exit(f'{info}\nPyOQP: dl-find meci optimization only support iopt 0–3, but found {iopt}') + + if ims not in [1, 2, 3]: + exit(f'{info}\nPyOQP: dl-find meci optimization only support ims 1-3, but found {ims}') + + if ims == 3 and icoord not in [10, 11, 12, 13, 14]: + exit(f'{info}\nPyOQP: dl-find meci optimization {ims} only support icoord 10–14, but found {icoord}') + + if ims in [1, 2] and icoord not in [0, 1, 2, 3, 4]: + exit(f'{info}\nPyOQP: dl-find meci optimization {ims} only support icoord 0–4, but found {icoord}') + + +def check_nac_input(config, info): + check_energy_input(config, info) + check_energy_input(config, info) + method = config['input']['method'] + td = config['tdhf']['type'] + nproc = config['nac']['nproc'] + ncpu = multiprocessing.cpu_count() + mpi = MPIManager().size + omp = int(os.environ['OMP_NUM_THREADS']) + + info += f'\n[input] method={method}' + + if method == 'hf': + exit(f'{info}\nPyOQP: scf cannot compute nac for hf calculations\n') + + if method == 'tdhf' and td != 'mrsf': + exit(f'{info}\nPyOQP: tdhf cannot compute nac for non-mrsf calculations\n') + + if MPIManager().use_mpi: + info += f'\nMPI process={mpi}' + info += f'\nOMP_NUM_THREADS={omp}' + print(f'{info}\nPyOQP: nac requested {mpi * omp} cpus! {ncpu} available on a single node') + else: + info += f'\n[nac] nproc={nproc}' + info += f'\nOMP_NUM_THREADS={omp}' + print(f'{info}\nPyOQP: nac requested {nproc * omp} cpus! {ncpu} available') + + if nproc < 1: + exit(f'{info}\nPyOQP: nac requested {nproc} process is illegal') + + +def check_soc_input(config, info): + not_available('soc', info) + + +def check_neb_input(config, info): + not_available('neb', info) + + +def check_hess_input(config, info): + check_energy_input(config, info) + method = config['input']['method'] + state = config['hess']['state'] + nproc = config['hess']['nproc'] + restart = config['hess']['restart'] + ncpu = multiprocessing.cpu_count() + mpi = MPIManager().size + omp = int(os.environ['OMP_NUM_THREADS']) + + info += f'\n[input] method={method}' + info += f'\n[hess] state={state}' + + if method == 'hf' and state > 0: + exit(f'{info}\nPyOQP: scf cannot compute hessian for state > 0\n') + + if method == 'tdhf' and state == 0: + exit(f'{info}\nPyOQP: tdhf cannot compute hessian for state = 0 \n') + + if restart != 'read': + if MPIManager().use_mpi: + info += f'\nMPI process={mpi}' + info += f'\nOMP_NUM_THREADS={omp}' + print(f'{info}\nPyOQP: hessian requested {mpi * omp} cpus! {ncpu} available on a single node') + else: + info += f'\n[hess] nproc={nproc}' + info += f'\nOMP_NUM_THREADS={omp}' + print(f'{info}\nPyOQP: hessian requested {nproc * omp} cpus! {ncpu} available') + + if nproc < 1: + exit(f'{info}\nPyOQP: hessian requested {nproc} process is illegal') + + +def check_nacme_input(config, info): + check_energy_input(config, info) + +def check_qmmm_input(config, info): + if not config['input']['qmmm_flag']: + exit(f'{info}\nPyOQP: QM/MM group reguires qmmm_flag=True in input group.') + diff --git a/pyoqp/oqp/utils/qmmm.py b/pyoqp/oqp/utils/qmmm.py new file mode 100644 index 0000000..46023ed --- /dev/null +++ b/pyoqp/oqp/utils/qmmm.py @@ -0,0 +1,802 @@ +import openmm as mm +import openmm.app as app +import openmm.unit as unit +import os +import numpy as np +import math + +# Information defined in OQP input +force_field = None +nonbondedMethod = None +constraints = None +rigidWater = False + +nSteps = 1 +timeStep = 1 +istate = 0 + +# List of QM atoms (defined in OQP) +qm_atoms = np.array((),dtype=np.uint64) + +# Initialize forcefield informations +pdb_file = None +pdb0 = None +forcefield0 = None +system0 = None + +# Initialize gradient information +gradient_qmmm = np.array((),dtype=np.float64) + +#Bonded info (generated here) +bonded_info = [] + +# For computing Link atom's g_factor: Define the array of covalent radii +covalent_radii = [ + 0.0, 0.354, 0.849, 1.336, 1.010, 0.838, + 0.757, 0.700, 0.658, 0.668, 0.920, 1.539, 1.421, + 1.244, 1.117, 1.101, 1.064, 1.044, 1.032, 1.953, + 1.761, 1.513, 1.412, 1.402, 1.345, 1.382, 1.270, + 1.241, 1.164, 1.302, 1.193, 1.260, 1.197, 1.211, + 1.190, 1.192, 1.147, 2.260, 2.052, 1.698, 1.564, + 1.473, 1.467, 1.322, 1.478, 1.332, 1.338, 1.386, + 1.403, 1.459, 1.398, 1.407, 1.386, 1.382, 1.267, + 2.570, 2.277, 1.943, 1.841, 1.823, 1.816, 1.801, + 1.780, 1.771, 1.735, 1.732, 1.710, 1.696, 1.673, + 1.660, 1.637, 1.671, 1.611, 1.511, 1.392, 1.372, + 1.372, 1.371, 1.364, 1.262, 1.340, 1.518, 1.459, + 1.512, 1.500, 1.545, 1.420, 2.880, 2.512, 1.983, + 1.721, 1.711, 1.684, 1.666, 1.657, 1.660, 1.801, + 1.761, 1.750, 1.724, 1.712, 1.689, 1.679, 1.698, + 1.850 +] + +#Conversion factor +kj_per_mol_to_hartree = 0.000380879 # 1 kJ/mol = 0.000380879 Hartree +bohr_to_nm = 0.052917721092 # 1 bohr = 0.052917721092 nm + +# Define unit charge and zero charge +unit_charge = 1.0 * unit.elementary_charge +zero_charge = 0.0 * unit.elementary_charge +zero_sigma = 0.0 * unit.nanometer +zero_epsilon = 0.0 * unit.kilojoule_per_mole + +#!---------------------------------------------------------------------- + +def _openmm_from_oqp(force_field,nonbondedMethod,constraints,water): + + if nonbondedMethod == 'NoCutoff': + electrostatics = app.NoCutoff + else: + exit(f"The {nonbondedMethod} type of electrostatics is not yet implemented") + + if constraints is None: + const = None + elif constraints == 'HBonds': + const = app.HBonds + elif constraints == 'AllBonds': + const = app.AllBonds + elif constraints == 'HAngles': + const = app.HAngles + else: + exit(f"Unknown type of constraint {constraints}, select a value between None, HBonds, AllBonds or HAngles") + + return force_field,electrostatics,const,water + +#!---------------------------------------------------------------------- + +def _openmm_atoms_bonded(topology, atom_index1, atom_index2): +# +# Purpose: Check if two atoms are bonded (needed for link atom placement) +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + for bond in topology.bonds(): + if (bond[0].index == atom_index1 and bond[1].index == atom_index2) or \ + (bond[0].index == atom_index2 and bond[1].index == atom_index1): + return True + return False + +#!---------------------------------------------------------------------- + +def _deactivate_bonded_all(system): +# +# Purpose: Deactivate bonded interactions between all atoms +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + #1) Deactivate all bonds + harmonic_bond_force = None + for force in system.getForces(): + if isinstance(force, mm.HarmonicBondForce): + harmonic_bond_force = force + break + if harmonic_bond_force is None: + raise ValueError("No HarmonicBondForce found in the system") + for i in range(harmonic_bond_force.getNumBonds()): + particle1, particle2, length, k = harmonic_bond_force.getBondParameters(i) + harmonic_bond_force.setBondParameters(i, particle1, particle2, length, 0.0 * unit.kilojoule_per_mole / unit.nanometer**2) + + #2) Deactivate all angles + harmonic_angle_force = None + for force in system.getForces(): + if isinstance(force, mm.HarmonicAngleForce): + harmonic_angle_force = force + break + if harmonic_angle_force is None: + raise ValueError("No HarmonicAngleForce found in the system") + for i in range(harmonic_angle_force.getNumAngles()): + particle1, particle2, particle3, angle, k = harmonic_angle_force.getAngleParameters(i) + harmonic_angle_force.setAngleParameters(i, particle1, particle2, particle3, angle, 0.0 * unit.kilojoule_per_mole / unit.radian**2) + + #3) Deactivate all torsions + periodic_torsion_force = None + for force in system.getForces(): + if isinstance(force, mm.PeriodicTorsionForce): + periodic_torsion_force = force + break + if periodic_torsion_force is None: + raise ValueError("No PeriodicTorsionForce found in the system") + for i in range(periodic_torsion_force.getNumTorsions()): + particle1, particle2, particle3, particle4, periodicity, phase, k = periodic_torsion_force.getTorsionParameters(i) + periodic_torsion_force.setTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, 0.0 * unit.kilojoule_per_mole) + +#!---------------------------------------------------------------------- + +def _deactivate_bonded_qm(system): +# +# Purpose: Deactivate bonded interactions between QM atoms +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + #1) Deactivate all QM-QM bonds + harmonic_bond_force = None + for force in system.getForces(): + if isinstance(force, mm.HarmonicBondForce): + harmonic_bond_force = force + break + if harmonic_bond_force is None: raise ValueError("No HarmonicBondForce found in the system") + + for i in range(harmonic_bond_force.getNumBonds()): + particle1, particle2, length, k = harmonic_bond_force.getBondParameters(i) + if particle1 in qm_atoms and particle2 in qm_atoms: + harmonic_bond_force.setBondParameters(i, particle1, particle2, length, 0.0 * unit.kilojoule_per_mole / unit.nanometer**2) + + #2) Deactivate all QM-QM-QM angles + harmonic_angle_force = None + for force in system.getForces(): + if isinstance(force, mm.HarmonicAngleForce): + harmonic_angle_force = force + break + if harmonic_angle_force is None: raise ValueError("No HarmonicAngleForce found in the system") + + for i in range(harmonic_angle_force.getNumAngles()): + particle1, particle2, particle3, angle, k = harmonic_angle_force.getAngleParameters(i) + if particle1 in qm_atoms and particle2 in qm_atoms and particle3 in qm_atoms: + harmonic_angle_force.setAngleParameters(i, particle1, particle2, particle3, angle, 0.0 * unit.kilojoule_per_mole / unit.radian**2) + + #3) Deactivate all QM-QM-QM-X or X-QM-QM-QM torsions (X=QM or MM) + periodic_torsion_force = None + for force in system.getForces(): + if isinstance(force, mm.PeriodicTorsionForce): + periodic_torsion_force = force + break + if periodic_torsion_force is None: raise ValueError("No PeriodicTorsionForce found in the system") + + for i in range(periodic_torsion_force.getNumTorsions()): + particle1, particle2, particle3, particle4, periodicity, phase, k = periodic_torsion_force.getTorsionParameters(i) + if (particle2 in qm_atoms and particle3 in qm_atoms) and (particle1 in qm_atoms or particle4 in qm_atoms): + periodic_torsion_force.setTorsionParameters(i, particle1, particle2, particle3, particle4, periodicity, phase, 0.0 * unit.kilojoule_per_mole) + +#!---------------------------------------------------------------------- + +def openmm_info(): +# +# Purpose: Update coordinates for OpenMM +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + + pdb = pdb0 + ff,electrostatics,const,water = _openmm_from_oqp(force_field,nonbondedMethod,constraints,rigidWater) + forcefield = app.ForceField(*ff) + system = forcefield.createSystem(pdb.topology, nonbondedMethod=electrostatics, constraints=const, rigidWater=water) + + xyz = [] + mass = [] + at_num = [] + natoms = system.getNumParticles() + +# Iterate over all particles and print their masses + for atom in pdb.topology.atoms(): + at_num.append(atom.element.atomic_number) + mass.append(system.getParticleMass(atom.index).value_in_unit(unit.dalton)) + pos=pdb.positions[atom.index].value_in_unit(unit.bohr) + xyz.append(pos.x) + xyz.append(pos.y) + xyz.append(pos.z) + + return xyz,mass,natoms,at_num + +#!---------------------------------------------------------------------- + +def openmm_dump(pdbfile=None): +# +# Purpose: Write PDB in filename +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + + app.PDBFile.writeFile(pdb0.topology, pdb0.positions, file=pdbfile) + for atom in pdb0.topology.atoms(): + i=atom.index + pos=pdb0.positions[i] + +#!---------------------------------------------------------------------- + +def openmm_update_system(coordinates): +# +# Purpose: Update coordinates for OpenMM +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + +############ +## Update coordinates +############ + xyz=coordinates.reshape((-1, 3)) + for atom in pdb0.topology.atoms(): + i=atom.index + pos=mm.Vec3(xyz[i][0],xyz[i][1],xyz[i][2])*unit.bohr + pdb0.positions[i]=pos + +#!---------------------------------------------------------------------- + +def openmm_init(atom_list): +# +# Purpose: Initialize atom list, forcefield, pdb and system +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + +### Check if atom list is ok + if len(atom_list) == 0: exit(f"\nError!! Atom list not defined!\n") + + pdb = pdb0 + +########### +## Read the PDB and load the forcefield information +############ + if pdb is None: +### The pdb file should be defined in the OQP input and it must exist + if not os.path.exists(pdb_file): exit(f"\nError!! PDB file is required in QM/MM runs!\n") + pdb = app.PDBFile(pdb_file) +# Attribute forcefield and create system from pdb topology + ff,electrostatics,const,water = _openmm_from_oqp(force_field,nonbondedMethod,constraints,rigidWater) + forcefield = app.ForceField(*ff) + system = forcefield.createSystem(pdb.topology, nonbondedMethod=electrostatics, constraints=const, rigidWater=water) + +#Shift atom list and verify it is not out of bounds + final_atom_list = np.array((),dtype=np.uint32) + for i in atom_list: + if int(i)-1 < 0 or int(i)-1 > system.getNumParticles(): + exit(f"\nError!! Atom list out of bounds!\n") + final_atom_list=np.append(final_atom_list,int(i)-1) + + return final_atom_list,pdb,forcefield,system + +#!---------------------------------------------------------------------- + +def openmm_system(): +# +# Purpose: Main function to pass all information OQP needs to know to perform QM/MM calculation +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# +### Reset data to be returned + num_atoms = 0 #Number of atoms (QM+link atoms) + q = [] #Atomic numbers (QM+link atoms) + mass = [] #Masses (QM+link atoms) in Daltons + x = [] #X positions (QM+link atoms) in Bohr + y = [] #Y positions (QM+link atoms) in Bohr + z = [] #Z positions (QM+link atoms) in Bohr + + pdb = pdb0 + forcefield = forcefield0 + system = system0 + +# Set up integrator + integrator = mm.LangevinIntegrator(300*unit.kelvin, 1/unit.picosecond, 0.002*unit.picoseconds) + simulation = app.Simulation(pdb.topology, system, integrator) + simulation.context.setPositions(pdb.positions) + +# Get the atom's position + positions = simulation.context.getState(getPositions=True).getPositions() +### Get the information about QM atoms (qm_info) and their connectivity (bonded_info, used to generate link atom list) + qm_info = [] + for qm_atom in pdb.topology.atoms(): + if qm_atom.index in qm_atoms: +# Retrieve atomic number, mass and position => save it in qm_info + qm_number = qm_atom.element.atomic_number + qm_mass = system.getParticleMass(qm_atom.index).value_in_unit(unit.dalton) + qm_position = positions[qm_atom.index].value_in_unit(unit.bohr) + qm_info.append((qm_number,qm_mass,[qm_position[0],qm_position[1],qm_position[2]])) +# Check if the specified atoms are bonded => save it in bonded_info + if not bonded_info: + for mm_atom in pdb.topology.atoms(): + if mm_atom.index not in qm_atoms: + if _openmm_atoms_bonded(pdb.topology, qm_atom.index, mm_atom.index): + mm_number = mm_atom.element.atomic_number + g_factor=(covalent_radii[1]+covalent_radii[qm_number])/(covalent_radii[qm_number]+covalent_radii[mm_number]) + bonded_info.append((qm_atom.index,mm_atom.index,g_factor)) + +### Compute total number of atoms (QM+Link Atoms) + num_atoms=len(qm_info)+len(bonded_info) + +### Adding qm_info to OQP's vectors + for i in qm_info: + q.append(i[0]) + mass.append(i[1]) + x.append(i[2][0]) + y.append(i[2][1]) + z.append(i[2][2]) + +### Adding Link atoms to OQP's vectors (Warning: For the moment, only hydrogen capping is allowed!) + linkatom_error = False + for i in range(len(bonded_info)): + if len(bonded_info[i]) != 0: + q.append(1) + mass.append(1.00782503223) + qm_index,mm_index,g_factor=bonded_info[i] + qm_position = positions[qm_index].value_in_unit(unit.bohr) + mm_position = positions[mm_index].value_in_unit(unit.bohr) + if g_factor >= 1.0 or g_factor <=0: + if not linkatom_error: + print("Error!! You should reconsider your QM/MM partitioning between:") + print(f" - QM({qm_atoms[i]}) and MM({mm_index})") + linkatom_error = True + x.append(qm_position[0]+g_factor*(mm_position[0]-qm_position[0])) + y.append(qm_position[1]+g_factor*(mm_position[1]-qm_position[1])) + z.append(qm_position[2]+g_factor*(mm_position[2]-qm_position[2])) + +# If error in link atom partitioning, exit! + if linkatom_error: exit() + + return num_atoms, x, y, z, q, mass + +#!---------------------------------------------------------------------- + +def openmm_potential(xyz): +# +# Purpose: Compute the MM potential on the xyz coordinates (QM+Link Atoms) +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + +### Reset data to be returned + mmpot = np.array(()) #MM potential (QM+link atoms) in Hartree/a.u. of charge + + pdb = pdb0 + ff,electrostatics,const,water = _openmm_from_oqp(force_field,nonbondedMethod,constraints,rigidWater) + forcefield = app.ForceField(*ff) + system = forcefield.createSystem(pdb.topology, nonbondedMethod=electrostatics, constraints=const, rigidWater=water) + +# forcefield = forcefield0 +# system = system0 + +#Focus on non-bonded interactions + forces = { force.__class__.__name__ : force for force in system.getForces() } + nonbonded_force = forces['NonbondedForce'] + if nonbonded_force is None: raise ValueError("No NonbondedForce found in the system") + +#Set all QM charges to zero and the VdW parameters to zero at first, including exceptions (scaling factors) + for index in range(nonbonded_force.getNumParticles()): + if index in qm_atoms: + nonbonded_force.setParticleParameters(index, zero_charge, zero_sigma, zero_epsilon) + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if p1 in qm_atoms or p2 in qm_atoms: + nonbonded_force.setExceptionParameters(i, p1, p2, zero_charge*zero_charge, zero_sigma, zero_epsilon) + +#Add link atoms to the list with a zero charge. Add exceptions inherited from QM bonded atom + num_real_particles=system.getNumParticles() + for ila in range(len(bonded_info)): + la_index=num_real_particles+ila + qm_index,mm_index,g_factor=bonded_info[ila] + system.addParticle(1.00782503223) + nonbonded_force.addParticle(zero_charge, zero_sigma, zero_epsilon) + link_atom_position=mm.Vec3(xyz[len(qm_atoms)+ila][0], + xyz[len(qm_atoms)+ila][1], + xyz[len(qm_atoms)+ila][2])*bohr_to_nm + pdb.positions.append(link_atom_position*unit.nanometer) + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if p1 == qm_index: + nonbonded_force.addException(la_index, p2, zero_charge*zero_charge, zero_sigma, zero_epsilon) + elif p2 == qm_index: + nonbonded_force.addException(p1, la_index, zero_charge*zero_charge, zero_sigma, zero_epsilon) + +# Create a simulation object with link atom + integrator = mm.LangevinIntegrator(300*unit.kelvin, 1/unit.picosecond, 0.002*unit.picoseconds) + simulation = app.Simulation(pdb.topology, system, integrator) + simulation.context.setPositions(pdb.positions) + +# Reset the pdb.positions to contain only real atoms + pdb.positions=pdb.positions[:num_real_particles] + +#Set positions and compute the total energy + simulation.context.reinitialize(preserveState=True) # Reinitialize to apply changes + state_without_charge = simulation.context.getState(getEnergy=True) + potential_energy_without_charge = state_without_charge.getPotentialEnergy() +# +# Potential on QM Atoms: Loop over all QM particles to get the potential at each QM center +# + for qm_index in qm_atoms: +#Get the potential with unit charge on target particle + nonbonded_force.setParticleParameters(qm_index, unit_charge, zero_sigma, zero_epsilon) + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if p1 == qm_index and p2 not in qm_atoms: + charge, sigma, epsilon = nonbonded_force.getParticleParameters(p2) + nonbonded_force.setExceptionParameters(i, p1, p2, charge*unit_charge, zero_sigma, zero_epsilon) + elif p2 == qm_index and p1 not in qm_atoms: + charge, sigma, epsilon = nonbonded_force.getParticleParameters(p1) + nonbonded_force.setExceptionParameters(i, p1, p2, charge*unit_charge, zero_sigma, zero_epsilon) + simulation.context.reinitialize(preserveState=True) # Reinitialize to apply changes + +# Compute the electrostatic potential energy contribution at the target particle + state_with_charge = simulation.context.getState(getEnergy=True) + potential_energy_with_charge = state_with_charge.getPotentialEnergy() + electrostatic_energy = potential_energy_with_charge - potential_energy_without_charge + electrostatic_energy_value = electrostatic_energy.value_in_unit(unit.kilojoule_per_mole) + +#Get the potential back to zero + nonbonded_force.setParticleParameters(qm_index, zero_charge, zero_sigma, zero_epsilon) + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if p1 == qm_index or p2 == qm_index: + nonbonded_force.setExceptionParameters(i, p1, p2, zero_charge*zero_charge, zero_sigma, zero_epsilon) + simulation.context.reinitialize(preserveState=True) # Reinitialize to apply changes + +# Print the potential at each QM center in Hartrees + mmpot=np.append(mmpot,[electrostatic_energy_value*kj_per_mol_to_hartree],axis=0) + +# +# Potential on Link Atoms: Loop over all link atoms to get the potential at each LA center +# + for ila in range(len(bonded_info)): + la_index = num_real_particles+ila + qm_index,mm_index,g_factor=bonded_info[ila] +#Get the potential with unit charge on target particle + nonbonded_force.setParticleParameters(la_index, unit_charge, zero_sigma, zero_epsilon) + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if p1 == la_index and p2 not in qm_atoms: + charge, sigma, epsilon = nonbonded_force.getParticleParameters(p2) + nonbonded_force.setExceptionParameters(i, la_index, p2, charge*unit_charge, zero_sigma, zero_epsilon) + elif p2 == la_index and p1 not in qm_atoms: + charge, sigma, epsilon = nonbonded_force.getParticleParameters(p1) + nonbonded_force.setExceptionParameters(i, p1, la_index, charge*unit_charge, zero_sigma, zero_epsilon) + simulation.context.reinitialize(preserveState=True) # Reinitialize to apply changes + +# Compute the electrostatic potential energy contribution at the target particle + state_with_charge = simulation.context.getState(getEnergy=True) + potential_energy_with_charge = state_with_charge.getPotentialEnergy() + electrostatic_energy = potential_energy_with_charge - potential_energy_without_charge + electrostatic_energy_value = electrostatic_energy.value_in_unit(unit.kilojoule_per_mole) + +#Get the potential back to zero + nonbonded_force.setParticleParameters(la_index, zero_charge, zero_sigma, zero_epsilon) + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if p1 == la_index and (p2 not in qm_atoms and p2 < num_real_particles): + nonbonded_force.setExceptionParameters(i, la_index, p2, zero_charge*zero_charge, zero_sigma, zero_epsilon) + elif p2 == la_index and (p1 not in qm_atoms and p1 < num_real_particles): + nonbonded_force.setExceptionParameters(i, p1, la_index, zero_charge*zero_charge, zero_sigma, zero_epsilon) + simulation.context.reinitialize(preserveState=True) # Reinitialize to apply changes + +# Print the potential at each QM center in Hartrees + mmpot=np.append(mmpot,[electrostatic_energy_value*kj_per_mol_to_hartree],axis=0) + + return mmpot + +#!---------------------------------------------------------------------- + +def openmm_energy(): +# +# Purpose: Get the MM energy components after deactivating QM-QM interactions +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + pdb = pdb0 + ff,electrostatics,const,water = _openmm_from_oqp(force_field,nonbondedMethod,constraints,rigidWater) + forcefield = app.ForceField(*ff) + system = forcefield.createSystem(pdb.topology, nonbondedMethod=electrostatics, constraints=const, rigidWater=water) + +# Deactivate all QM-QM bonded interactions and QM bonded terms including exceptions + _deactivate_bonded_qm(system) + forces = { force.__class__.__name__ : force for force in system.getForces() } + nonbonded_force = forces['NonbondedForce'] + if nonbonded_force is None: raise ValueError("No NonbondedForce found in the system") + + for i in range(nonbonded_force.getNumParticles()): + charge, sigma, epsilon = nonbonded_force.getParticleParameters(i) + if i in qm_atoms: + nonbonded_force.setParticleParameters(i, charge*0.0, sigma, epsilon) + + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if p1 in qm_atoms or p2 in qm_atoms: + nonbonded_force.setExceptionParameters(i, p1, p2, chargeProd*0.0, sigma*0.0, epsilon*0.0) + +# Enable energy decomposition + for i, f in enumerate(system.getForces()): + f.setForceGroup(i) + +# Set up integrator + integrator = mm.LangevinIntegrator(300*unit.kelvin, 1/unit.picosecond, 0.002*unit.picoseconds) + simulation = app.Simulation(pdb.topology, system, integrator) + simulation.context.setPositions(pdb.positions) + +# Energy decomposition + energy=np.zeros((24)) + for i, f in enumerate(system.getForces()): + state = simulation.context.getState(getEnergy=True, groups={i}) + current_energy=state.getPotentialEnergy().value_in_unit(unit.kilojoule_per_mole)*kj_per_mol_to_hartree + key=f.getName() +# +# This code only works for python version >= 3.10 +# +# match key: +# case "HarmonicBondForce": energy[0]=current_energy +# case "HarmonicAngleForce": energy[1]=current_energy +# case "PeriodicTorsionForce": energy[2]=current_energy +# case "RBTorsionForce": energy[3]=current_energy +# case "NonbondedForce": energy[4]=current_energy +# case "CustomBondForce": energy[5]=current_energy +# case "CustomAngleForce": energy[6]=current_energy +# case "CustomTorsionForce": energy[7]=current_energy +# case "CustomNonbondedForce": energy[8]=current_energy +# case "CustomExternalForce": energy[9]=current_energy +# case "GBSAOBCForce": energy[10]=current_energy +# case "CustomGBForce": energy[11]=current_energy +# case "AmoebaMultipoleForce": energy[12]=current_energy +# case "AmoebaVdwForce": energy[13]=current_energy +# case "AmoebaBondForce": energy[14]=current_energy +# case "AmoebaAngleForce": energy[15]=current_energy +# case "AmoebaTorsionForce": energy[16]=current_energy +# case "AmoebaOutOfPlaneBendForce": energy[17]=current_energy +# case "AmoebaPiTorsionForce": energy[18]=current_energy +# case "AmoebaStretchBendForce": energy[19]=current_energy +# case "AmoebaUreyBradleyForce": energy[20]=current_energy +# case "AmoebaVdw14Force": energy[21]=current_energy +# case "CMMotionRemover": energy[22]=current_energy +# case "CustomCompoundBondForce": energy[23]=current_energy +# case "MonteCarloBarostat": energy[24]=current_energy +# case _: +# exit(f"Unknown {key} Force. Please update the qmmm.py code!") +# +# Substituted by if/elif equivalent: + if key == "HarmonicBondForce": energy[0] = current_energy + elif key == "HarmonicAngleForce": energy[1] = current_energy + elif key == "PeriodicTorsionForce": energy[2] = current_energy + elif key == "RBTorsionForce": energy[3] = current_energy + elif key == "NonbondedForce": energy[4] = current_energy + elif key == "CustomBondForce": energy[5] = current_energy + elif key == "CustomAngleForce": energy[6] = current_energy + elif key == "CustomTorsionForce": energy[7] = current_energy + elif key == "CustomNonbondedForce": energy[8] = current_energy + elif key == "CustomExternalForce": energy[9] = current_energy + elif key == "GBSAOBCForce": energy[10] = current_energy + elif key == "CustomGBForce": energy[11] = current_energy + elif key == "AmoebaMultipoleForce": energy[12] = current_energy + elif key == "AmoebaVdwForce": energy[13] = current_energy + elif key == "AmoebaBondForce": energy[14] = current_energy + elif key == "AmoebaAngleForce": energy[15] = current_energy + elif key == "AmoebaTorsionForce": energy[16] = current_energy + elif key == "AmoebaOutOfPlaneBendForce": energy[17] = current_energy + elif key == "AmoebaPiTorsionForce": energy[18] = current_energy + elif key == "AmoebaStretchBendForce": energy[19] = current_energy + elif key == "AmoebaUreyBradleyForce": energy[20] = current_energy + elif key == "AmoebaVdw14Force": energy[21] = current_energy + elif key == "CMMotionRemover": energy[22] = current_energy + elif key == "CustomCompoundBondForce": energy[23] = current_energy + elif key == "MonteCarloBarostat": energy[24] = current_energy + else: + exit(f"Unknown {key} Force. Please update the qmmm.py code!") + +# Remove QM-QM electrostatics +# Deactivate all bonded interactions + _deactivate_bonded_all(system) + +# Deactivate all QM-MM and MM-MM nonbonded interactions + for index in range(nonbonded_force.getNumParticles()): + if index not in qm_atoms: + nonbonded_force.setParticleParameters(index, 0.0, 0.0, 0.0) + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if (p1 not in qm_atoms and p2 not in qm_atoms) or (p1 in qm_atoms and p2 not in qm_atoms) or (p1 not in qm_atoms and p2 in qm_atoms): + nonbonded_force.setExceptionParameters(i, p1, p2, chargeProd*0.0, sigma*0.0, epsilon*0.0) + simulation.context.reinitialize(preserveState=True) # Reinitialize to apply changes + state = simulation.context.getState(getEnergy=True) + for i, f in enumerate(system.getForces()): + state = simulation.context.getState(getEnergy=True, groups={i}) + current_energy=state.getPotentialEnergy().value_in_unit(unit.kilojoule_per_mole)*kj_per_mol_to_hartree + key=f.getName() +# +# This code only works for python version >= 3.10 +# +# match key: +# case "HarmonicBondForce": energy[0]-=current_energy +# case "HarmonicAngleForce": energy[1]-=current_energy +# case "PeriodicTorsionForce": energy[2]-=current_energy +# case "RBTorsionForce": energy[3]-=current_energy +# case "NonbondedForce": energy[4]-=current_energy +# case "CustomBondForce": energy[5]-=current_energy +# case "CustomAngleForce": energy[6]-=current_energy +# case "CustomTorsionForce": energy[7]-=current_energy +# case "CustomNonbondedForce": energy[8]-=current_energy +# case "CustomExternalForce": energy[9]-=current_energy +# case "GBSAOBCForce": energy[10]-=current_energy +# case "CustomGBForce": energy[11]-=current_energy +# case "AmoebaMultipoleForce": energy[12]-=current_energy +# case "AmoebaVdwForce": energy[13]-=current_energy +# case "AmoebaBondForce": energy[14]-=current_energy +# case "AmoebaAngleForce": energy[15]-=current_energy +# case "AmoebaTorsionForce": energy[16]-=current_energy +# case "AmoebaOutOfPlaneBendForce": energy[17]-=current_energy +# case "AmoebaPiTorsionForce": energy[18]-=current_energy +# case "AmoebaStretchBendForce": energy[19]-=current_energy +# case "AmoebaUreyBradleyForce": energy[20]-=current_energy +# case "AmoebaVdw14Force": energy[21]-=current_energy +# case "CMMotionRemover": energy[22]-=current_energy +# case "CustomCompoundBondForce": energy[23]-=current_energy +# case "MonteCarloBarostat": energy[24]-=current_energy +# case _: +# exit(f"Unknown {key} Force. Please update the qmmm.py code!") +# +# Substituted by if/elif equivalent: + if key == "HarmonicBondForce": energy[0] -= current_energy + elif key == "HarmonicAngleForce": energy[1] -= current_energy + elif key == "PeriodicTorsionForce": energy[2] -= current_energy + elif key == "RBTorsionForce": energy[3] -= current_energy + elif key == "NonbondedForce": energy[4] -= current_energy + elif key == "CustomBondForce": energy[5] -= current_energy + elif key == "CustomAngleForce": energy[6] -= current_energy + elif key == "CustomTorsionForce": energy[7] -= current_energy + elif key == "CustomNonbondedForce": energy[8] -= current_energy + elif key == "CustomExternalForce": energy[9] -= current_energy + elif key == "GBSAOBCForce": energy[10] -= current_energy + elif key == "CustomGBForce": energy[11] -= current_energy + elif key == "AmoebaMultipoleForce": energy[12] -= current_energy + elif key == "AmoebaVdwForce": energy[13] -= current_energy + elif key == "AmoebaBondForce": energy[14] -= current_energy + elif key == "AmoebaAngleForce": energy[15] -= current_energy + elif key == "AmoebaTorsionForce": energy[16] -= current_energy + elif key == "AmoebaOutOfPlaneBendForce": energy[17] -= current_energy + elif key == "AmoebaPiTorsionForce": energy[18] -= current_energy + elif key == "AmoebaStretchBendForce": energy[19] -= current_energy + elif key == "AmoebaUreyBradleyForce": energy[20] -= current_energy + elif key == "AmoebaVdw14Force": energy[21] -= current_energy + elif key == "CMMotionRemover": energy[22] -= current_energy + elif key == "CustomCompoundBondForce": energy[23] -= current_energy + elif key == "MonteCarloBarostat": energy[24] -= current_energy + else: + exit(f"Unknown {key} Force. Please update the qmmm.py code!") + + return energy + +#!---------------------------------------------------------------------- + +def openmm_gradient(xyz,oqp_charges): +# +# Purpose: Get the MM gradient components after deactivating QM-QM interactions +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + pdb = pdb0 + ff,electrostatics,const,water = _openmm_from_oqp(force_field,nonbondedMethod,constraints,rigidWater) + forcefield = app.ForceField(*ff) + system = forcefield.createSystem(pdb.topology, nonbondedMethod=electrostatics, constraints=const, rigidWater=water) + +# Remove QM-QM bonded interactions + _deactivate_bonded_qm(system) + +# Add QM charges to the atom list + forces = { force.__class__.__name__ : force for force in system.getForces() } + nonbonded_force = forces['NonbondedForce'] + if nonbonded_force is None: raise ValueError("No NonbondedForce found in the system") + for index in range(nonbonded_force.getNumParticles()): + if index in qm_atoms: + charge, sigma, epsilon = nonbonded_force.getParticleParameters(index) + charge = oqp_charges[np.where(qm_atoms == index)[0][0]]*unit.elementary_charge + nonbonded_force.setParticleParameters(index, charge, sigma, epsilon) + +#Add link atoms to the list with a zero charge. Add exceptions inherited from QM bonded atom + num_real_particles=system.getNumParticles() + for ila in range(len(bonded_info)): + la_index=num_real_particles+ila + qm_index,mm_index,g_factor=bonded_info[ila] + system.addParticle(1.00782503223) + nonbonded_force.addParticle(oqp_charges[len(qm_atoms)+ila],0.0,0.0) + link_atom_position=mm.Vec3(xyz[len(qm_atoms)+ila][0], + xyz[len(qm_atoms)+ila][1], + xyz[len(qm_atoms)+ila][2])*bohr_to_nm + pdb.positions.append(link_atom_position*unit.nanometer) + +# Create a simulation object with link atom + integrator = mm.LangevinIntegrator(300*unit.kelvin, 1/unit.picosecond, 0.002*unit.picoseconds) + simulation = app.Simulation(pdb.topology, system, integrator) + simulation.context.setPositions(pdb.positions) + +# Reset the pdb.positions to contain only real atoms + pdb.positions=pdb.positions[:num_real_particles] + +# Compute gradient + state = simulation.context.getState(getForces=True) + gradient = -1.0*state.getForces(asNumpy=True).value_in_unit(unit.kilojoule_per_mole/unit.bohr)*kj_per_mol_to_hartree + +# Remove QM-QM non-bonded + +# Deactivate all bonded interactions + _deactivate_bonded_all(system) + +# Deactivate all QM-MM and MM-MM nonbonded interactions, including exceptions + for index in range(nonbonded_force.getNumParticles()): + if index not in qm_atoms: + nonbonded_force.setParticleParameters(index, 0.0, 0.0, 0.0) + + for i in range(nonbonded_force.getNumExceptions()): + p1, p2, chargeProd, sigma, epsilon = nonbonded_force.getExceptionParameters(i) + if (p1 not in qm_atoms and p2 not in qm_atoms) or (p1 in qm_atoms and p2 not in qm_atoms) or (p1 not in qm_atoms and p2 in qm_atoms): + nonbonded_force.setExceptionParameters(i, p1, p2, chargeProd*0.0, sigma*0.0, epsilon*0.0) + + simulation.context.reinitialize(preserveState=True) # Reinitialize to apply changes + state = simulation.context.getState(getForces=True) + gradient2 = -1.0*state.getForces(asNumpy=True).value_in_unit(unit.kilojoule_per_mole/unit.bohr)*kj_per_mol_to_hartree + + gradient_mm = gradient - gradient2 + + grad_qm=np.zeros((len(qm_atoms)+len(bonded_info),3)) + + for i in qm_atoms: + grad_qm[np.where(i == qm_atoms)[0][0],:]=gradient_mm[i,:] + + for ila in range(len(bonded_info)): + grad_qm[len(qm_atoms)+ila,:]=gradient_mm[num_real_particles+ila,:] + + return grad_qm,gradient_mm + +#!---------------------------------------------------------------------- + +def form_gradient_qmmm(gradient_qm,gradient_mm): +# +# Purpose: Form the final QM/MM gradient +# Author: Miquel Huix-Rotllant +# Date: 31/07/2024 +# + pdb = pdb0 + ff,electrostatics,const,water = _openmm_from_oqp(force_field,nonbondedMethod,constraints,rigidWater) + forcefield = app.ForceField(*ff) + system = forcefield.createSystem(pdb.topology, nonbondedMethod=electrostatics, constraints=const, rigidWater=water) + +# For link atoms, first project back the contributions + num_real_particles=system.getNumParticles() + +# Project back link atom contributions to real atoms + if len(bonded_info) != 0: + for i in range(len(bonded_info)): + qm_index,mm_index,g_factor=bonded_info[i] + gradient_qm[0,np.where(qm_index == qm_atoms)[0][0],:]+=(1.0-g_factor)*gradient_qm[0,len(qm_atoms)+i,:] + gradient_mm[mm_index,:]+=g_factor*gradient_qm[0,len(qm_atoms)+i,:] + +# Form the final QM/MM gradient + qmmm_gradient=np.zeros((num_real_particles,3),dtype=np.float64) + for i in range(num_real_particles): + if i not in qm_atoms: + qmmm_gradient[i,:]=gradient_mm[i,:] + elif i in qm_atoms: + qmmm_gradient[i,:]=gradient_qm[0,np.where(i == qm_atoms)[0][0],:] + else: + exit(f"Error in form_gradient_qmmm, unknown atom index {i}") + + return qmmm_gradient diff --git a/pyoqp/setup.py b/pyoqp/setup.py index 3f24e9e..f37d34d 100644 --- a/pyoqp/setup.py +++ b/pyoqp/setup.py @@ -19,7 +19,8 @@ 'libdlfind>=0.0.3', # 'dftd4>=3.5.0', 'cffi>=1.16.0', - 'mpi4py>=4.0.0', +# 'mpi4py>=4.0.0', + 'openmm>=8.1.1', ], extras_require={ }, diff --git a/source/integrals/grd1.F90 b/source/integrals/grd1.F90 index 9fedd3c..b7113c4 100644 --- a/source/integrals/grd1.F90 +++ b/source/integrals/grd1.F90 @@ -31,6 +31,7 @@ module grd1 public grad_ee_kinetic public grad_en_hellman_feynman public grad_en_pulay + public grad_elpot contains @@ -385,7 +386,122 @@ SUBROUTINE grad_ee_kinetic(basis, denab, de, logtol) DEALLOCATE(de_priv) END SUBROUTINE +!MHR START +!------------------------------------------------------------------------------- + +!> @brief Electrostatic potential grid derivative contributions +!> @details Compute derivative integrals of type +!> @note No relativistic methods available +! +!> @author Vladimir Mironov, Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Apr, 2024_ QMMM modifications +!> +!> @param[in] coord grid coordinates +!> @param[in] zq grid weights +!> @param[in,out] denab density matrix in packed format, remains unchanged on return +!> @param[in] l2 dimension of density matrix array +!> @param[out] de integrals + SUBROUTINE grad_elpot(basis, coord, zq, denab, de, logtol) + + REAL(kind=dp), INTENT(INOUT) :: denab(:) + type(basis_set), intent(inout) :: basis + real(kind=dp), contiguous, intent(in) :: coord(:) + real(kind=dp), intent(in) :: zq + + REAL(kind=dp) :: de(:,:) + + INTEGER :: l2, & + ii, jj + + REAL(kind=dp), optional :: logtol + + REAL(kind=dp) :: de1(3) + REAL(kind=dp), ALLOCATABLE :: de_priv(:,:), dens(:,:) + + REAL(kind=dp) :: dernuc(3), tol + LOGICAL :: out, dbg, norm + + TYPE(shell_t) :: shi, shj + TYPE(shpair_t) :: cntp + + INTEGER :: nat + + dbg = .false. + out = .false. + + IF (dbg) WRITE(iw,'(/10X,38(1H-)/10X,"GRADIENT INCLUDING AO DERIVATIVE TERMS"/10X,38(1H-))') + + nat = ubound(de, 2) + l2 = basis%nbf + + !IF (dbg) write(iw,*) "OMP 1E GRD (TVDER)", basis%nshell, nat + + if (present(logtol)) then + tol = logtol + else + tol = tol_default + end if + + norm = .true. + allocate(dens(basis%nbf,basis%nbf), source=0.0d0) + call unpack_matrix(denab, dens) + IF (norm) THEN + CALL bas_norm_matrix(dens, basis%bfnrm, basis%nbf) + END IF + +! temporary storage for 1e gradient + ALLOCATE(de_priv, mold=de) + de_priv = 0.0d0 + +!$omp parallel & +!$omp private( & +!$omp ii, jj, & +!$omp shi, shj, cntp, & +!$omp dernuc, & +!$omp de1 & +!$omp ) & +!$omp reduction(+:de_priv) + + CALL cntp%alloc(basis) + +!$omp do schedule(dynamic) +! I shell + DO ii = 1, basis%nshell + + + CALL shi%fetch_by_id(basis, ii) + de1 = 0.0 + +! J shell + DO jj = 1, basis%nshell + + CALL shj%fetch_by_id(basis, jj) + + CALL cntp%shell_pair(basis, shi, shj, tol, dup=.false.) + IF (cntp%numpairs==0) CYCLE + +! Nuclear attraction derivative + CALL comp_coulomb_der1(cntp, coord(:), -zq, dens(basis%ao_offset(ii):, basis%ao_offset(jj):), dernuc) + de1 = de1 + 2*dernuc(1:3) +! End of primitive loops + + END DO + + de_priv(:,shi%atid) = de_priv(:,shi%atid) + de1(:) + + END DO +!$omp end do +! End of shell loops +!$omp end parallel + + de = de + de_priv + DEALLOCATE(de_priv) + + END SUBROUTINE +!MHR END !------------------------------------------------------------------------------- !> @brief Basis function derivative contributions to gradient diff --git a/source/integrals/int1.F90 b/source/integrals/int1.F90 index 9bfb1e3..102e3bf 100644 --- a/source/integrals/int1.F90 +++ b/source/integrals/int1.F90 @@ -45,6 +45,7 @@ module int1 private public omp_hst + public omp_qmmm public multipole_integrals public electrostatic_potential public basis_overlap @@ -177,6 +178,59 @@ subroutine omp_hst(basis, coord, zq, h, s, t, z, debug, logtol) !------------------------------------------------------------------------------- +!> @brief Driver for conventional ESP QM/MM integrals on a grid around QM atoms +! +!> @details Compute one electron integrals and core Hamiltonian, +!> - V is evaluated by Gauss-Rys quadrature, then \f$ h = T+V \f$ +!> Also, do \f$ L_z \f$ integrals for atoms and linear cases. +!> Also, do FMO ESP integrals if needed. +!> This subroutine is capable to do integrals in parallel using +!> both OpenMP and MPI. It it helpful when running large FMO jobs. +! +!> @note Based on `HSANDT` subroutine from file `INT1.SRC` +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Jul, 2024_ Initial release +!> +!> @param[in] i QM center index +!> @param[in] ttt ESP integral weight +!> @param[in,out] chg_op one-electron ESP atomic charge operator in packet format + subroutine omp_qmmm(basis, i, coord, ttt, chg_op, nat, logtol) + + use precision, only: dp + use basis_tools, only: basis_set + + type(basis_set), intent(in) :: basis + real(real64), contiguous, intent(in) :: coord(:,:), ttt(:,:) + real(real64), contiguous, intent(inout) :: chg_op(:) + real(real64), optional, intent(in) :: logtol + integer, intent(in) :: i, nat + real(real64) :: tol + + integer :: l1, l2 + + if (present(logtol)) then + tol = logtol + else + tol = log(10.0_dp)*20 + end if + +! Exclude 1e potential in ESDIM, because density is used, +! not point charges + + l1 = basis%nbf + l2 = l1*(l1+1)/2 + + chg_op(:)=0 + call nuc_ints(basis, coord, ttt(i,:), chg_op(:), tol) + call bas_norm_matrix(chg_op(:), basis%bfnrm, l1) + + end subroutine + +!------------------------------------------------------------------------------- + !> @brief Driver for multipole integrals ! !> @details Compute one electron multipole integrals @@ -726,7 +780,7 @@ subroutine nuc_ints(basis, coord, zq, h, tol) call cntp%alloc(basis) -! I shell +! I shell do ii = basis%nshell, 1, -1 call shi%fetch_by_id(basis, ii) @@ -1015,6 +1069,39 @@ SUBROUTINE int1_kin_ovl(cntp, dokinetic, sblk, tblk) !-------------------------------------------------------------------------------- +!> @brief Compute contracted block of Coulomb 1e integrals +!> @param[in] cntp shell pair data +!> @param[in] xyz coordinates of particles +!> @param[in] c charges of particles +!> @param[in] nat number of particles +!> @param[in] chgtol cut-off for charge +!> @param[inout] blk block of 1e Coulomb integrals +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Jul, 2024_ Initial release +! + SUBROUTINE int1_coul_xyz_u(cntp, xyz, c, blk) +!dir$ attributes inline :: int1_coulxyz_u + TYPE(shpair_t), INTENT(IN) :: cntp + REAL(REAL64), CONTIGUOUS, INTENT(IN) :: xyz(:) + REAL(REAL64), INTENT(IN) :: c + REAL(REAL64), CONTIGUOUS, INTENT(INOUT) :: blk(:) + + INTEGER :: ig + +!dir$ assume_aligned blk : 64 + +! Interaction with point charge + DO ig = 1, cntp%numpairs + CALL comp_coulomb_int1_prim(cntp, ig, xyz, -c, blk) + END DO + + END SUBROUTINE + +!-------------------------------------------------------------------------------- + !> @brief Compute contracted block of Coulomb 1e integrals !> @param[in] cntp shell pair data !> @param[in] xyz coordinates of particles diff --git a/source/mathlib/lapack_wrap.F90 b/source/mathlib/lapack_wrap.F90 index 0febba9..75e699c 100644 --- a/source/mathlib/lapack_wrap.F90 +++ b/source/mathlib/lapack_wrap.F90 @@ -85,6 +85,92 @@ subroutine oqp_dgesv_i64(n, nrhs, a, lda, ipiv, b, ldb, info) end subroutine +!MHR START +!---------------------------------------------------------------------- +!> Compute the inverse of a real matrix using LAPACK routine dgetri. +!! @param n Size of the matrix. +!! @param a On entry, the matrix to be inverted. On exit, the inverted +!matrix. +!! @param lda Leading dimension of the array `a`. +!! @param ipiv Integer array of size `n` containing pivot indices computed by +!dgetrf. +!! @param work Workspace array of size `lwork`. +!! @param lwork Size of the workspace array. +!! @param info On exit, info = 0 for successful exit. If info = -i, the i-th +!argument had an illegal value. +subroutine oqp_dgetri_i64(n, a, lda, ipiv, work, lwork, info) + integer :: n, lda, ipiv(:), lwork, info + real(dp) :: a(lda,*), work(lwork) + integer(blas_int) :: n_, lda_, lwork_, info_ + logical :: ok + ! Check arguments + if (ARG_CHECK) then + ok = .true. + ! Check if the size of n, lda, and lwork is within acceptable bounds + ok = ok .and. (n <= HUGE_BLAS_INT) + ok = ok .and. (lda <= HUGE_BLAS_INT) + ok = ok .and. (lwork <= HUGE_BLAS_INT) + ! If any of the arguments are invalid, display error message and abort + if (.not. ok) call show_message(ERRMSG, WITH_ABORT) + end if + ! Convert integer arguments to blas_int type + n_ = int(n, blas_int) + lda_ = int(lda, blas_int) + lwork_ = int(lwork, blas_int) + ! Call LAPACK routine dgetri to compute the inverse of the matrix + call dgetri(n_, a, lda_, ipiv, work, lwork_, info_) + ! Assign the output info_ to the info argument + info = info_ +end subroutine +!> LU factorization of a general matrix using LAPACK routine dgetrf. +!! This subroutine computes the LU factorization of a general M-by-N matrix A +!using partial pivoting with row interchanges. +!! The factorization has the form: +!! A = P * L * U +!! where P is a permutation matrix, L is lower triangular with unit diagonal +!elements (lower trapezoidal if m > n), +!! and U is upper triangular (upper trapezoidal if m < n). +!! This subroutine overwrites the input matrix A with its factors L and U. +!! @param m Number of rows in the matrix A. +!! @param n Number of columns in the matrix A. +!! @param a On entry, the matrix to be factorized. On exit, the factors L +!and U. +!! @param lda Leading dimension of the array `a`, must be at least max(1,m). +!! @param ipiv Integer array of dimension at least min(m,n) containing the +!pivot indices. +!! @param info On exit, info = 0 for successful exit. If info = -i, the i-th +!argument had an illegal value. +subroutine oqp_dgetrf_i64(m, n, a, lda, ipiv, info) + integer :: m, n, lda, ipiv(min(m,n)), info + real(dp) :: a(lda,*) + integer(blas_int) :: m_, n_, lda_, info_ + integer(blas_int), allocatable :: ipiv_(:) + logical :: ok + ! Check arguments + if (ARG_CHECK) then + ok = .true. + ! Check if the size of m, n, and lda is within acceptable bounds + ok = ok .and. (m <= HUGE_BLAS_INT) + ok = ok .and. (n <= HUGE_BLAS_INT) + ok = ok .and. (lda <= HUGE_BLAS_INT) + ! If any of the arguments are invalid, display error message and abort + if (.not. ok) call show_message(ERRMSG, WITH_ABORT) + end if + ! Convert integer arguments to blas_int type + m_ = int(m, blas_int) + n_ = int(n, blas_int) + lda_ = int(lda, blas_int) + ! Allocate memory for 'ipiv_' and copy 'ipiv' to it + allocate(ipiv_(min(m_,n_))) + ipiv_ = int(ipiv, blas_int) + ! Call LAPACK routine dgetrf to perform LU factorization + call dgetrf(m_, n_, a, lda_, ipiv_, info_) + ! Assign the output info_ to the info argument + info = info_ + ! Free allocated memory +end subroutine oqp_dgetrf_i64 +!MHR END + !---------------------------------------------------------------------- subroutine oqp_dsysv_i64(uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info) @@ -291,46 +377,4 @@ subroutine oqp_dtrttp_i64(uplo, n, a, lda, ap, info) !---------------------------------------------------------------------- -!> Compute the inverse of a real matrix using LAPACK routine dgetri. -!! @param n Size of the matrix. -!! @param a On entry, the matrix to be inverted. On exit, the inverted matrix. -!! @param lda Leading dimension of the array `a`. -!! @param ipiv Integer array of size `n` containing pivot indices computed by dgetrf. -!! @param work Workspace array of size `lwork`. -!! @param lwork Size of the workspace array. -!! @param info On exit, info = 0 for successful exit. If info = -i, the i-th argument had an illegal value. -subroutine oqp_dgetri_i64(n, a, lda, ipiv, work, lwork, info) - integer :: n, lda, ipiv(:), lwork, info - real(dp) :: a(lda,*), work(lwork) - - integer(blas_int) :: n_, lda_, lwork_, info_ - logical :: ok - - ! Check arguments - if (ARG_CHECK) then - ok = .true. - - ! Check if the size of n, lda, and lwork is within acceptable bounds - ok = ok .and. (n <= HUGE_BLAS_INT) - ok = ok .and. (lda <= HUGE_BLAS_INT) - ok = ok .and. (lwork <= HUGE_BLAS_INT) - - ! If any of the arguments are invalid, display error message and abort - if (.not. ok) call show_message(ERRMSG, WITH_ABORT) - end if - - ! Convert integer arguments to blas_int type - n_ = int(n, blas_int) - lda_ = int(lda, blas_int) - lwork_ = int(lwork, blas_int) - - ! Call LAPACK routine dgetri to compute the inverse of the matrix - call dgetri(n_, a, lda_, ipiv, work, lwork_, info_) - - ! Assign the output info_ to the info argument - info = info_ -end subroutine - -!---------------------------------------------------------------------- - end module diff --git a/source/modules/hf_gradient.F90 b/source/modules/hf_gradient.F90 index dd58753..0a3e835 100644 --- a/source/modules/hf_gradient.F90 +++ b/source/modules/hf_gradient.F90 @@ -149,6 +149,7 @@ subroutine hf_1e_grad(infos, basis) use constants, only: tol_int use grd1, only: eijden, grad_nn, grad_ee_overlap, & grad_ee_kinetic, grad_en_hellman_feynman, grad_en_pulay + use qmmm_mod, only: grad_esp_qmmm implicit none @@ -211,6 +212,9 @@ subroutine hf_1e_grad(infos, basis) ! Pulay force call grad_en_pulay(basis, xyz, zn, dens, grad, logtol=tol) +! QM/MM force + if(infos%control%qmmm_flag) call grad_esp_qmmm(infos, dens, grad,logtol=tol) + end associate end subroutine diff --git a/source/modules/int1e.F90 b/source/modules/int1e.F90 index 51bc09c..27e06a3 100644 --- a/source/modules/int1e.F90 +++ b/source/modules/int1e.F90 @@ -33,6 +33,7 @@ subroutine int1e(infos) use strings, only: Cstring, fstring use physical_constants, only: BOHR_TO_ANGSTROM use printing, only: print_module_info + use qmmm_mod, only: oqp_esp_qmmm implicit none @@ -46,9 +47,11 @@ subroutine int1e(infos) ! tagarray real(kind=dp), contiguous, pointer :: & - hcore(:), tmat(:), smat(:) + hcore(:), tmat(:), smat(:), Hqmmm(:), mm_potential(:) character(len=*), parameter :: tags_general(3) = (/ character(len=80) :: & - OQP_SM, OQP_TM, OQP_Hcore /) + OQP_SM, OQP_TM, OQP_Hcore/) + character(len=*), parameter :: tags_qmmm(2) = (/ character(len=80) :: & + OQP_mm_potential, OQP_Hqmmm/) logical dbg dbg = .false. @@ -78,13 +81,13 @@ subroutine int1e(infos) ! Allocate H, S and T matrices nbf2 = basis%nbf*(basis%nbf+1)/2 + nat = ubound(infos%atoms%zn,1) +! If you remove records, you removes the tagarray of OQP_mm_potential defined in Python. call infos%dat%remove_records(tags_general) - call infos%dat%reserve_data(OQP_SM, TA_TYPE_REAL64, nbf2, comment=OQP_SM_comment) call infos%dat%reserve_data(OQP_TM, TA_TYPE_REAL64, nbf2, comment=OQP_TM_comment) call infos%dat%reserve_data(OQP_Hcore, TA_TYPE_REAL64, nbf2, comment=OQP_Hcore_comment) - call data_has_tags(infos%dat, tags_general, module_name, subroutine_name, WITH_ABORT) call tagarray_get_data(infos%dat, OQP_SM, smat) call tagarray_get_data(infos%dat, OQP_TM, tmat) @@ -98,8 +101,31 @@ subroutine int1e(infos) tol = log(10.0d0)*tol_int call omp_hst(basis, infos%atoms%xyz, infos%atoms%zn, hcore, smat, tmat, logtol=tol) +! Compute QM/MM interaction + if(infos%control%qmmm_flag) then + call infos%dat%reserve_data(OQP_Hqmmm, TA_TYPE_REAL64, nbf2, comment=OQP_Hqmmm_comment) + call data_has_tags(infos%dat, tags_qmmm, module_name, subroutine_name, WITH_ABORT) + call tagarray_get_data(infos%dat, OQP_Hqmmm, Hqmmm) + call tagarray_get_data(infos%dat, OQP_mm_potential, mm_potential) + + write(iw,"(/1X,' Computing ESP One Electron Integrals (QM/MM) '/)") + write(iw,"('External MM potential:'/)") + do i=1,nat + write(iw,"(i4,1X,f12.8)") i, mm_potential(i) + end do +! Compute QM/MM contribution to core Hamiltonian + call oqp_esp_qmmm(infos, Hqmmm, mm_potential, smat, logtol=tol) +! Add QM/MM contribution to core Hamiltonian + hcore = hcore + hqmmm + write(iw,"(/1X,' ... End of ESP One Electron Integrals ... '/)") + endif + if (dbg) then - write(iw,'(/"BARE NUCLEUS HAMILTONIAN INTEGRALS (H=T+V)")') + if(infos%control%qmmm_flag) then + write(iw,'(/"BARE NUCLEUS HAMILTONIAN INTEGRALS (H=T+V+QM/MM)")') + else + write(iw,'(/"BARE NUCLEUS HAMILTONIAN INTEGRALS (H=T+V)")') + end if call print_sym_labeled(hcore,nbf, basis) write(iw,'(/"OVERLAP MATRIX")') @@ -107,6 +133,11 @@ subroutine int1e(infos) write(iw,'(/"KINETIC ENERGY INTEGRALS")') call print_sym_labeled(tmat,nbf, basis) + + if (infos%control%qmmm_flag) then + write(iw,'(/"QM/MM HAMILTONIAN INTEGRALS")') + call print_sym_labeled(Hqmmm,nbf, basis) + end if end if write(iw,"(/1X,'...... End Of One Electron Integrals ......'/)") diff --git a/source/modules/qmmm.F90 b/source/modules/qmmm.F90 new file mode 100644 index 0000000..b02e668 --- /dev/null +++ b/source/modules/qmmm.F90 @@ -0,0 +1,710 @@ +module qmmm_mod + + use oqp_linalg + use resp_mod, only: add_atom_grid + implicit none + + character(len=*), parameter :: module_name = "qmmm_mod" + + private + public get_mm_energy + public form_esp_charges + public print_mm_energy + public oqp_esp_qmmm + public grad_esp_qmmm + + contains + +!> @brief Compute the MM and atomic QM/MM contributions to energy +! +!> @detail Classical contributions to energy and the QM atomic charge times MM potential +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Oct, 2024_ Initial release + subroutine get_mm_energy(infos,emm) + use precision, only: dp + use oqp_tagarray_driver + use types, only: information + use messages, only: WITH_ABORT + implicit none + + character(len=*), parameter :: subroutine_name = "get_mm_energy" + + real(kind=dp), intent(inout) :: emm + real(kind=dp), contiguous, pointer :: mm_potential(:), mm_energy(:) + character(len=*), parameter :: tags_qmmm(2) = (/ character(len=80) :: & + OQP_mm_potential, OQP_mm_energy /) + + type(information), target, intent(inout) :: infos + + emm = 0.0d0 + if(.not.infos%control%qmmm_flag) return + + call data_has_tags(infos%dat, tags_qmmm, module_name, subroutine_name, WITH_ABORT) + call tagarray_get_data(infos%dat, OQP_mm_energy, mm_energy) + call tagarray_get_data(infos%dat, OQP_mm_potential, mm_potential) + emm=sum(mm_energy)+dot_product(mm_potential,infos%atoms%zn) + + return + + end subroutine get_mm_energy + +!-------------------------------------------------------------------------------- + +!> @brief Print MM energy decomposition in file +! +!> @detail Classical contributions to energy and the QM atomic charge times MM potential +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Oct, 2024_ Initial release + subroutine print_mm_energy(infos) + use precision, only: dp + use oqp_tagarray_driver + use types, only: information + use messages, only: WITH_ABORT + use io_constants, only: iw + implicit none + + character(len=*), parameter :: subroutine_name = "get_mm_energy" + + real(kind=dp), contiguous, pointer :: mm_potential(:), mm_energy(:),partial_charges(:) + character(len=*), parameter :: tags_qmmm(3) = (/ character(len=80) :: & + OQP_mm_potential, OQP_mm_energy, OQP_partial_charges /) + + type(information), target, intent(inout) :: infos + + if(.not.infos%control%qmmm_flag) return + + call data_has_tags(infos%dat, tags_qmmm, module_name, subroutine_name, WITH_ABORT) + call tagarray_get_data(infos%dat, OQP_mm_energy, mm_energy) + call tagarray_get_data(infos%dat, OQP_mm_potential, mm_potential) + call tagarray_get_data(infos%dat, OQP_partial_charges, partial_charges) + + write(iw,*) +! QM/MM interaction energy + write(iw,"(' QM/MM: Electrostatic energy = ',F20.10)")& + dot_product(partial_charges,mm_potential) +! Classical forcefields + if(abs(mm_energy(1)).gt.1.0d-10)& + write(iw,"(' MM: Bonded energy = ',F20.10)") mm_energy(1) + if(abs(mm_energy(2)).gt.1.0d-10)& + write(iw,"(' MM: Angle bending energy = ',F20.10)") mm_energy(2) + if(abs(mm_energy(3)).gt.1.0d-10)& + write(iw,"(' MM: Periodic Torsion energy = ',F20.10)") mm_energy(3) + if(abs(mm_energy(4)).gt.1.0d-10)& + write(iw,"(' MM: RBTorsion energy = ',F20.10)") mm_energy(4) + if(abs(mm_energy(5)).gt.1.0d-10)& + write(iw,"(' MM: Nonbonded energy = ',F20.10)") mm_energy(5) + +! Custom Classical forcefields + if(abs(mm_energy(6)).gt.1.0d-10)& + write(iw,"(' MM: CustomBonded energy = ',F20.10)") mm_energy(6) + if(abs(mm_energy(7)).gt.1.0d-10)& + write(iw,"(' MM: CustomAngle energy = ',F20.10)") mm_energy(7) + if(abs(mm_energy(8)).gt.1.0d-10)& + write(iw,"(' MM: CustomTorsion energy = ',F20.10)") mm_energy(8) + if(abs(mm_energy(9)).gt.1.0d-10)& + write(iw,"(' MM: CustomNonbonded energy = ',F20.10)") mm_energy(9) + if(abs(mm_energy(10)).gt.1.0d-10)& + write(iw,"(' MM: CustomExternal energy = ',F20.10)") mm_energy(10) +! Implicit solvation + if(abs(mm_energy(11)).gt.1.0d-10)& + write(iw,"(' MM: GBSAOBC energy = ',F20.10)") mm_energy(11) + if(abs(mm_energy(12)).gt.1.0d-10)& + write(iw,"(' MM: CustomGB energy = ',F20.10)") mm_energy(12) + +! Amoeba forcefields + if(abs(mm_energy(13)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaMultipole energy = ',F20.10)") mm_energy(13) + if(abs(mm_energy(14)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaVdw energy = ',F20.10)") mm_energy(14) + if(abs(mm_energy(15)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaBond energy = ',F20.10)") mm_energy(15) + if(abs(mm_energy(16)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaAngle energy = ',F20.10)") mm_energy(16) + if(abs(mm_energy(17)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaTorsion energy = ',F20.10)") mm_energy(17) + if(abs(mm_energy(18)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaOutOfPlaneBend energy = ',F20.10)") mm_energy(18) + if(abs(mm_energy(19)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaPiTorsion energy = ',F20.10)") mm_energy(19) + if(abs(mm_energy(20)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaStretchBend energy = ',F20.10)") mm_energy(20) + if(abs(mm_energy(21)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaUreyBradley energy = ',F20.10)") mm_energy(21) + if(abs(mm_energy(22)).gt.1.0d-10)& + write(iw,"(' MM: AmoebaVdw14 energy = ',F20.10)") mm_energy(22) +! Extras + if(abs(mm_energy(23)).gt.1.0d-10)& + write(iw,"(' MM: CMMotion energy = ',F20.10)") mm_energy(23) + if(abs(mm_energy(24)).gt.1.0d-10)& + write(iw,"(' MM: CustomCompoundBond energy = ',F20.10)") mm_energy(24) + if(abs(mm_energy(25)).gt.1.0d-10)& + write(iw,"(' MM: MonteCarloBarostad = ',F20.10)") mm_energy(25) + write(iw,*) + + call print_charges(infos,partial_charges,iw) + + return + + end subroutine print_mm_energy + +!-------------------------------------------------------------------------------- + +!> @brief Compute the MM and atomic QM/MM contributions to energy +! +!> @detail Classical contributions to energy and the QM atomic charge times MM potential +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Oct, 2024_ Initial release + + subroutine form_esp_charges(infos,dmat,nbf) + use precision, only: dp + use oqp_tagarray_driver + use basis_tools, only: basis_set + use messages, only: show_message, with_abort + use types, only: information + use strings, only: Cstring, fstring + use int1, only: omp_qmmm + use constants, only: tol_int + use mathlib, only: traceprod_sym_packed + implicit none + + character(len=*), parameter :: subroutine_name = "form_esp_charges" + + integer, intent(in) :: nbf + real(kind=dp), intent(in) :: dmat(nbf*nbf) + + type(information), target, intent(inout) :: infos + + real(kind=dp), contiguous, pointer :: partial_charges(:) + character(len=*), parameter :: tags_qmmm(1) = (/ character(len=80) :: & + OQP_partial_charges /) + + integer :: nat,nelec,nbf2,ok + type(basis_set), pointer :: basis + real(kind=dp), allocatable :: chg_op(:) + real(kind=dp), target, allocatable :: ttt(:,:), xyz(:,:) + real(kind=dp), pointer :: coord(:,:) + real(kind=dp), pointer :: wn(:,:) + + integer :: npt, nptcur + integer :: i + real(kind=dp) :: chg,corr + + integer, parameter :: nlayers = 4 + real(kind=dp), parameter :: & + layers(nlayers) = [1.4, 1.6, 1.8, 2.0] + integer, parameter :: npt_layer(nlayers) = [132,152,192,350] + integer, parameter :: typ_layer(nlayers) = [3,3,3,0] + + real(kind=dp) :: tol + + if(.not.infos%control%qmmm_flag) return + +! Load basis set + basis => infos%basis + basis%atoms => infos%atoms + +! Allocate memory + nbf2 = nbf*(nbf+1)/2 + nat = ubound(infos%atoms%zn,1) + nelec = infos%mol_prop%nelec + npt = nat * sum(npt_layer) + + allocate(xyz(npt,3), & + ttt(nat,npt), & + chg_op(nbf2), & + stat=ok) + if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT) + + call form_espf_grid(nat,npt,nlayers,layers,npt_layer,typ_layer,infos%atoms%zn,infos%atoms%xyz,xyz,ttt,nptcur) + +! Compute integrals and form partial charges + wn => ttt(:,:nptcur) + tol = log(10.0d0)*tol_int + + call infos%dat%reserve_data(OQP_partial_charges, TA_TYPE_REAL64, infos%mol_prop%natom, comment=OQP_partial_charges_comment) + call data_has_tags(infos%dat, tags_qmmm, module_name, subroutine_name, WITH_ABORT) + call tagarray_get_data(infos%dat, OQP_partial_charges, partial_charges) + +! Form partial charges and compute the correction for total charge + corr=nelec/nat + do i=1,nat + call omp_qmmm(basis, i, transpose(xyz(1:nptcur,:)), wn, chg_op, nat, logtol=tol) + chg = traceprod_sym_packed(dmat,chg_op,nbf) + partial_charges(i) = infos%atoms%zn(i) + chg + corr = corr + chg/nat + end do + +! Correct for the total charge + do i=1,nat + partial_charges(i) = partial_charges(i) - corr + end do + + return + + end subroutine form_esp_charges + +!-------------------------------------------------------------------------------- + +!> @brief Compute the QM/MM one-electron hamiltonian using ESPF method +!> @param[in] infos OQP handle +!> @param[in,out] hqmmm QM/MM one-electron hamiltonian +!> @param[in] mm_potential classical MM potential on QM centers +!> @param[in] smat overlap matrix +!> @param[in] logtol tolerance threshold for integrals +! +!> @detail This subroutine computes the one-electron hamiltonian for introducing the QM/MM interaction +!> in the core hamiltonian using the electrostatic potential fitted (ESPF) method. +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Jul, 2024_ Initial release + subroutine oqp_esp_qmmm(infos, Hqmmm, mm_potential, smat, logtol) + use precision, only: dp + use basis_tools, only: basis_set + use messages, only: show_message, with_abort + use types, only: information + use strings, only: Cstring, fstring + use int1, only: omp_qmmm + use lebedev, only: lebedev_get_grid + use elements, only: ELEMENTS_VDW_RADII + implicit none + + character(len=*), parameter :: subroutine_name = "oqp_esp_qmmm" + + type(information), target, intent(inout) :: infos + real(kind=dp), contiguous, intent(inout) :: Hqmmm(:) + real(kind=dp), contiguous, intent(in) :: smat(:),mm_potential(:) + real(kind=dp), intent(in) :: logtol + + integer :: nbf, nbf2, ok + type(basis_set), pointer :: basis + real(kind=dp), allocatable :: wt(:), chg_op(:) + real(kind=dp), allocatable :: xyz(:,:) + real(kind=dp), target, allocatable :: ttt(:,:) + real(kind=dp), pointer :: coord(:,:) + real(kind=dp), pointer :: wn(:,:) + + integer :: nat, npt, nptcur + integer :: i + real(kind=dp) :: mm_pot_av + + integer, parameter :: nlayers = 4 + real(kind=dp), parameter :: & + layers(nlayers) = [1.4, 1.6, 1.8, 2.0] + integer, parameter :: npt_layer(nlayers) = [132,152,192,350] + integer, parameter :: typ_layer(nlayers) = [3,3,3,0] + + logical :: restr + +! Load basis set + basis => infos%basis + basis%atoms => infos%atoms + +! Allocate memory + nbf = basis%nbf + nbf2 = nbf*(nbf+1)/2 + nat = ubound(infos%atoms%zn, 1) + npt = nat * sum(npt_layer) + + allocate(xyz(npt,3), & + wt(npt), & + ttt(nat,npt), & + chg_op(nbf2), & + stat=ok) + if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT) + +! Compute the numerical Lebedev grid (xyz) and integral weights (TTT) + call form_espf_grid(nat,npt,nlayers,layers,npt_layer,typ_layer,infos%atoms%zn,infos%atoms%xyz,xyz,ttt,nptcur) + +! Compute integrals and form QM/MM hamiltonian + wn => ttt(:,:nptcur) + mm_pot_av = sum(mm_potential)/nat + hqmmm = -smat*mm_pot_av + do i=1,nat + call omp_qmmm(basis, i, transpose(xyz(1:nptcur,:)), wn, chg_op, nat, logtol) + hqmmm = hqmmm + (mm_potential(i)-mm_pot_av)*chg_op + end do + + end subroutine oqp_esp_qmmm + +!-------------------------------------------------------------------------------- + +!> @brief Compute the gradient contributions from the QM/MM interaction energy using ESPF method +!> @param[in] infos OQP handle +!> @param[in] dens density matrix +!> @param[in,out] grad energy gradient +!> @param[in] logtol tolerance threshold for integrals +! +!> @detail This subroutine computes the analytic derivatives of the QM/MM interaction energy +!> from the electrostatic potential fitted (ESPF) method. Here the polarizable term +!> (Q^x*mm_potential) and the classical term (Q*mm_potential^x) computed by OpenMM +!> are added to the QM gradient. The MM atoms gradient is added in the python layer. +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Jul, 2024_ Initial release + subroutine grad_esp_qmmm(infos, dens, grad, logtol) + use oqp_tagarray_driver + use precision, only: dp + use basis_tools, only: basis_set + use messages, only: show_message, with_abort + use types, only: information + use strings, only: Cstring, fstring + use grd1, only: grad_elpot, grad_ee_overlap + use lebedev, only: lebedev_get_grid + use elements, only: ELEMENTS_VDW_RADII + implicit none + + character(len=*), parameter :: subroutine_name = "grad_esp_qmmm" + + type(information), target, intent(inout) :: infos + real(kind=dp), intent(inout) :: grad(:,:) + real(kind=dp), intent(inout) :: dens(:) + real(kind=dp), intent(in) :: logtol + + integer :: nbf, nbf2, ok + type(basis_set), pointer :: basis + real(kind=dp), allocatable :: wt(:) + real(kind=dp), allocatable :: xyz(:,:) + real(kind=dp), target, allocatable :: ttt(:,:) + + integer :: nat, npt, nptcur + integer :: i, j, k + + integer, parameter :: nlayers = 4 + real(kind=dp), parameter :: & + layers(nlayers) = [1.4, 1.6, 1.8, 2.0] + integer, parameter :: npt_layer(nlayers) = [132,152,192,350] + integer, parameter :: typ_layer(nlayers) = [3,3,3,0] + +!tagarray + real(kind=dp), contiguous, pointer :: partial_charges(:), mm_potential(:),& + gradient_mm(:,:) + character(len=*), parameter :: tags_qmmm(3) = (/ character(len=80) :: & + OQP_partial_charges, OQP_mm_potential, OQP_mm_gradient/) + + real(kind=dp) :: mm_pot_av + +! Load basis set + basis => infos%basis + basis%atoms => infos%atoms + +! Allocate memory + nbf = basis%nbf + nbf2 = nbf*(nbf+1)/2 + nat = ubound(infos%atoms%zn, 1) + npt = nat * sum(npt_layer) + + allocate(xyz(npt,3), & + wt(npt), & + ttt(nat,npt), & + stat=ok) + if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT) + + call form_espf_grid(nat,npt,nlayers,layers,npt_layer,typ_layer,infos%atoms%zn,infos%atoms%xyz,xyz,ttt,nptcur) + +! ESP gradient contribution +! Tagarray + call data_has_tags(infos%dat, tags_qmmm, module_name, subroutine_name, WITH_ABORT) + call tagarray_get_data(infos%dat, OQP_mm_potential, mm_potential) + call tagarray_get_data(infos%dat, OQP_partial_charges, partial_charges) + call tagarray_get_data(infos%dat, OQP_mm_gradient, gradient_mm) + +! Compute integrals and form ESP operators + +! Compute the corrected mm potential + mm_pot_av = sum(mm_potential)/nat + do i=1,nat + mm_potential(i)=mm_pot_av-mm_potential(i) + end do + +! Add integral gradient term, mm_potential*[(T^+T)^-1*T^+]*V^x + do i=1,nptcur + wt(i)=-dot_product(ttt(:,i),mm_potential) + call grad_elpot(basis, xyz(i,:), wt(i), dens, grad) + end do + +! Add overlap derivative correction for the total charge conservation + if(abs(mm_pot_av).gt.1.0e-6) then + dens=-mm_pot_av*dens + call grad_ee_overlap(basis, dens, grad, logtol) + dens=-dens/mm_pot_av + end if + +! Add weights gradient term, -mm_potential*[(T^+T)^-1*T^+]*T^xQ + q*mm_potential^x + call espf_grad(& + x=xyz(:nptcur,1),& + y=xyz(:nptcur,2),& + z=xyz(:nptcur,3),& + at=infos%atoms%xyz,& + wt=wt,& + zn=infos%atoms%zn,& + pchg=partial_charges,& + grad_mm=gradient_mm,& + grad=grad) + +! Restablish mm potential to the original one + do i=1,nat + mm_potential(i)=-mm_pot_av-mm_potential(i) + end do + + end subroutine grad_esp_qmmm + +!-------------------------------------------------------------------------------- + +!> @brief Add to the gradient the integral weight derivatives and classical MM contributions +!> @param[in] infos OQP handle +!> @param[in] dens density matrix +!> @param[in,out] grad energy gradient +!> @param[in] logtol tolerance threshold for integrals +! +!> @detail This subroutine computes the analytic derivatives of the QM/MM interaction energy +!> corresponding to the integral weight derivatives and the classical MM contributions +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Jul, 2024_ Initial release + subroutine espf_grad(x, y, z, at, wt, zn, pchg, grad_mm, grad) + use precision, only: dp + implicit none + + real(kind=dp), intent(in) :: x(:), y(:), z(:), at(:,:), zn(:) + real(kind=dp), intent(inout) :: grad(:,:) + real(kind=dp), contiguous, intent(in) :: pchg(:), grad_mm(:,:) + real(kind=dp), intent(inout) :: wt(:) + real(kind=dp), target, allocatable :: gradient_mm(:,:) + + integer :: i, j, nat, npts, dim1,dim2 + + nat = ubound(at, 2) + npts = ubound(x, 1) + + dim1=ubound(grad_mm,1) + dim2=ubound(grad_mm,2) + allocate(gradient_mm(dim2,dim1)) + gradient_mm=reshape(grad_mm, (/dim2,dim1/)) + + ! Compute the matrix of the inverse distances + do i = 1, nat + grad(1,i) = grad(1,i) + gradient_mm(1,i) + grad(2,i) = grad(2,i) + gradient_mm(2,i) + grad(3,i) = grad(3,i) + gradient_mm(3,i) + do j = 1, npts + grad(1,i) = grad(1,i) & + + wt(j)*(pchg(i)-zn(i))*(at(1,i)-x(j))/norm2(at(:,i)-[x(j),y(j),z(j)])**3 + grad(2,i) = grad(2,i) & + + wt(j)*(pchg(i)-zn(i))*(at(2,i)-y(j))/norm2(at(:,i)-[x(j),y(j),z(j)])**3 + grad(3,i) = grad(3,i) & + + wt(j)*(pchg(i)-zn(i))*(at(3,i)-z(j))/norm2(at(:,i)-[x(j),y(j),z(j)])**3 + end do + end do + + end subroutine espf_grad + +!-------------------------------------------------------------------------------- + +!> @brief Computes the integral weights for computing the ESPF integrals +!> @param[in] x,y,z Grid Cartesian coordinates +!> @param[in,out] ttt Integral weights +!> @param[in] at Atom Cartesian coordinates +! +!> @detail Form the electrostatic kernel T(nat,npts) and forms the pseudoinverse +!> (T^+T)^-1*T^+ (in which + = dagger) +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Jul, 2024_ Initial release + subroutine espf_weights(x, y, z, ttt, at) + use io_constants, only: iw + use precision, only: dp + implicit none + + real(kind=dp), intent(in) :: x(:), y(:), z(:), at(:,:) + real(kind=dp), intent(inout) :: ttt(:,:) + + integer :: i, j, nat, npts, lwork, info + real(kind=dp), allocatable :: b(:,:), r(:,:) + real(kind=dp), allocatable :: work(:) + real(kind=dp) :: workk(1) + + nat = ubound(at, 2) + npts = ubound(x, 1) + + allocate(b(npts,npts), r(nat,npts)) + + ! Compute the matrix of the inverse distances + b=0.0d0 + do j = 1, npts + b(j,j)=1.0d0 + do i = 1, nat + r(i,j) = 1/norm2(at(:,i) - [x(j), y(j), z(j)]) + end do + end do + + ! Workspace query for optimal lwork size + lwork = -1 + call dgels('T', nat, npts, npts, r, nat, b, npts, workk, lwork, info) + lwork = int(workk(1)) + allocate(work(lwork)) + + ! Solve the least-squares problem T * X ≈ B + call dgels('T', nat, npts, npts, r, nat, b, npts, work, lwork, info) + + ttt(:nat,:npts)=b(:nat,:npts) + + end subroutine espf_weights + +!-------------------------------------------------------------------------------- + +!> @brief Adds atomic-centered spherical grid to form the molecular grid for ESP calculations +!> @param[in] nat Number of QM centers +!> @param[in] npt Maximum number of gridpoints +!> @param[in] nlayers Number of layers for the grid +!> @param[in] layers Radius for the different layers +!> @param[in] npt_layer Number of points per layer +!> @param[in] typ_layer +!> @param[in] zn atomic charges +!> @param[in] atoms_xyz atomic coordinates +!> @param[in,out] xyz grid coordinates +!> @param[in,out] ttt grid weights +!> @param[in,out] nptcur final number of grid points +! +!> @detail Form the numerical grid for ESPF computations and construct the integral weights +!> This routine is adapted from oqp/modules/resp.F90 +! +!> @author Miquel Huix-Rotllant +! +! REVISION HISTORY: +!> @date _Oct, 2024_ Initial release + subroutine form_espf_grid(nat,npt,nlayers,layers,npt_layer,typ_layer,zn,atoms_xyz,xyz,ttt,nptcur) + use precision, only: dp + use lebedev, only: lebedev_get_grid + use elements, only: ELEMENTS_VDW_RADII + use messages, only: show_message, with_abort + implicit none + + integer, intent(in) :: nat,npt,nlayers + real(kind=dp), intent(in) :: zn(nat),atoms_xyz(3,nat),layers(nlayers) + integer, intent(in) :: npt_layer(nlayers),typ_layer(nlayers) + real(kind=dp), intent(inout) :: xyz(npt,3),ttt(nat,npt) + integer, intent(inout) :: nptcur + + real(kind=dp), allocatable :: leb(:,:), lebw(:) + real(kind=dp), allocatable :: vdwrad(:), wt(:) + integer, allocatable :: neigh(:) + + integer :: nadd, nleb, ok + integer :: i, j, k, layer + + allocate(leb(maxval(npt_layer),3), & + lebw(maxval(npt_layer)), & + vdwrad(nat), & + neigh(nat), & + wt(npt), & + stat=ok) + if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT) + +! Set up the grid + nptcur = 0 ! Current number of point in a grid + +! Loop over layers and add spherical grid points on each atom to the molecular grid + do layer = 1, nlayers + +! Get the atomic radii on which to place new grid layer + vdwrad = ELEMENTS_VDW_RADII(int(zn))*layers(layer) + +! Get grid + nleb = npt_layer(layer) + leb = 0 + lebw = 0 + call lebedev_get_grid(nleb, leb, lebw, typ_layer(layer)) + +! Add new grid layer for each atom, remove inner points + do i = 1, nat + call add_atom_grid( & + x=xyz(nptcur+1:,1), & + y=xyz(nptcur+1:,2), & + z=xyz(nptcur+1:,3), & + wts=wt(nptcur+1:), & + nadd=nadd, & + atpts=leb(:nleb,:), & + atwts=lebw(:nleb), & + atoms_xyz=atoms_xyz, & + atoms_rad=vdwrad, & + cur_atom=i, & + neighbours=neigh) + nptcur = nptcur + nadd + end do + end do + +! Compute the weights + call espf_weights( & + x=xyz(:nptcur,1), & + y=xyz(:nptcur,2), & + z=xyz(:nptcur,3), & + ttt=ttt, & + at = atoms_xyz) + + end subroutine form_espf_grid + +!-------------------------------------------------------------------------------- + +!> @brief Print partial ESPF charges +!> +!> @param[in] infos OQP handle +!> @param[in] chg atomic partial charges +! +!> @author Miquel Huix-Rotllant +! +!> @detail This is a modified copy of print_charges of oqp/modules/resp.F90 +! +! REVISION HISTORY: +!> @date _Oct, 2024_ Initial release + subroutine print_charges(infos,partial_charges,iw) + use precision, only: dp + use oqp_tagarray_driver + use elements, only: ELEMENTS_SHORT_NAME + use types, only: information + use messages, only: show_message, with_abort + implicit none + + type(information), intent(in) :: infos + real(kind=dp), intent(in), contiguous :: partial_charges(:) + integer, intent(in) :: iw + + integer :: i, elem, nat + + nat = ubound(infos%atoms%zn, 1) + + write(iw,'(2/)') + write(iw,'(4x,a)') '==============================' + write(iw,'(4x,a)') 'ESPF QM/MM charges calculation' + write(iw,'(4x,a)') '==============================' + + write(iw,'(/,30("^"))') + write(iw,'(/a8,a8,a14)') '#', 'Name', 'Charge' + write(iw,'(30("-"))') + + do i = 1, nat + elem = nint(infos%atoms%zn(i)) + write(iw,'(i8,a8,f14.6)') i, ELEMENTS_SHORT_NAME(elem), partial_charges(i) + end do + + write(iw,'(30("="))') + + end subroutine print_charges + +end module qmmm_mod diff --git a/source/modules/resp.F90 b/source/modules/resp.F90 index dc5612a..f92ebaf 100644 --- a/source/modules/resp.F90 +++ b/source/modules/resp.F90 @@ -8,6 +8,7 @@ module resp_mod private public oqp_resp_charges + public add_atom_grid public resp_charges_C !-------------------------------------------------------------------------------- @@ -223,7 +224,7 @@ subroutine oqp_resp_charges(infos) q=infos%atoms%zn, & pot=pot(:nptcur)) -! Fit ESP c8harges +! Fit ESP charges call chg_fit_mk( & x=xyz(:nptcur,1), & y=xyz(:nptcur,2), & diff --git a/source/scf.F90 b/source/scf.F90 index a1e0904..0446521 100644 --- a/source/scf.F90 +++ b/source/scf.F90 @@ -31,6 +31,7 @@ subroutine scf_driver(basis, infos, molGrid) use basis_tools, only: basis_set use scf_converger, only: scf_conv_result, scf_conv, & conv_cdiis, conv_ediis + use qmmm_mod, only: get_mm_energy,form_esp_charges,print_mm_energy implicit none @@ -41,7 +42,7 @@ subroutine scf_driver(basis, infos, molGrid) type(dft_grid_t), intent(in) :: molGrid integer :: i, ii, iter, nschwz, nbf, nbf_tri, nbf2, ok, maxit real(kind=dp) :: ehf, ehf1, nenergy, etot, diffe, e_old, psinrm, & - scalefactor,vne, vnn, vee, vtot, virial, tkin + scalefactor,vne, vnn, vee, vtot, virial, tkin, emm real(kind=dp), allocatable :: tempvec(:), lwrk(:), lwrk2(:) real(kind=dp), allocatable, target :: smat_full(:,:), pdmat(:,:), pfock(:,:), rohf_bak(:) real(kind=dp), allocatable, target :: dold(:,:), fold(:,:) @@ -171,6 +172,9 @@ subroutine scf_driver(basis, infos, molGrid) call tagarray_get_data(infos%dat, OQP_VEC_MO_B, mo_b) endif +! Compute the (fixed) classical MM contribution to energy if QM/MM run (otherwise emm=0.0d0) + call get_mm_energy(infos,emm) + allocate(smat_full(nbf,nbf), pdmat(nbf_tri,nfocks), pfock(nbf_tri,nfocks), & qmat(nbf,nbf), & stat=ok, & @@ -379,6 +383,9 @@ subroutine scf_driver(basis, infos, molGrid) etot = etot + eexc end if +!Adding MM energy + if(infos%control%qmmm_flag) etot = etot + emm + ! Forming ROHF Fock by combing Alpha and Beta Focks. if (scf_type == scf_rohf) then rohf_bak = pfock(:,1) @@ -545,6 +552,15 @@ subroutine scf_driver(basis, infos, molGrid) mo_energy_b = mo_energy_a end select + ! Construct ESPF partial charges and print MM energy in output (only done if QM/MM run) + select case (scf_type) + case (scf_rhf) + call form_esp_charges(infos,dmat_a,nbf) + case (scf_uhf,scf_rohf) + call form_esp_charges(infos,dmat_a+dmat_b,nbf) + end select + call print_mm_energy(infos) + ! Print out the molecular orbitals call print_mo_range(basis, infos, mostart=1, moend=nbf) diff --git a/source/tagarray_driver.F90 b/source/tagarray_driver.F90 index 836d931..6a0dc21 100644 --- a/source/tagarray_driver.F90 +++ b/source/tagarray_driver.F90 @@ -45,6 +45,11 @@ module oqp_tagarray_driver character(len=*), parameter, public :: OQP_nac = OQP_prefix // "nac" character(len=*), parameter, public :: OQP_td_states_phase = OQP_prefix // "td_states_phase" character(len=*), parameter, public :: OQP_td_states_overlap = OQP_prefix // "td_states_overlap" + character(len=*), parameter, public :: OQP_mm_potential = OQP_prefix // "mm_potential" + character(len=*), parameter, public :: OQP_Hqmmm = OQP_prefix // "hamiltonian_qmmm" + character(len=*), parameter, public :: OQP_partial_charges = OQP_prefix // "partial_charges" + character(len=*), parameter, public :: OQP_mm_energy = OQP_prefix // "mm_energy" + character(len=*), parameter, public :: OQP_mm_gradient = OQP_prefix // "mm_gradient" character(len=*), parameter, public :: OQP_DM_A_comment = "Alpha-spin triangle Density matrix" character(len=*), parameter, public :: OQP_DM_B_comment = "Beta-spin triangle Density matrix" @@ -66,6 +71,11 @@ module oqp_tagarray_driver character(len=*), parameter, public :: OQP_td_xpy_comment = OQP_prefix // "(X+Y) vector for target state in TD-DFT calculations" character(len=*), parameter, public :: OQP_td_xmy_comment = OQP_prefix // "(X-Y) vector for target state in TD-DFT calculations" character(len=*), parameter, public :: OQP_td_energies_comment = OQP_prefix // "Responce energies" + character(len=*), parameter, public :: OQP_mm_potential_comment = "MM potential" + character(len=*), parameter, public :: OQP_partial_charges_comment = "QM partial charges" + character(len=*), parameter, public :: OQP_mm_energy_comment = "MM energy" + character(len=*), parameter, public :: OQP_mm_gradient_comment = "MM gradient" + character(len=*), parameter, public :: OQP_Hqmmm_comment = "triangle QM/MM Hamiltonian matrix" character(len=*), parameter, public :: OQP_log_filename_comment = OQP_prefix // "log filename" character(len=*), parameter, public :: OQP_basis_filename_comment = OQP_prefix // "basis filename" character(len=*), parameter, public :: OQP_hbasis_filename_comment = OQP_prefix // "Huckel basis_filename for Huckel Guess" @@ -75,14 +85,15 @@ module oqp_tagarray_driver character(len=*), parameter, public :: OQP_td_states_phase_comment = OQP_prefix // "Bvecs phase sign with respect to Bvec_old" character(len=*), parameter, public :: OQP_td_states_overlap_comment = OQP_prefix // "Bvecs phase sign with respect to Bvec_old" character(len=*), parameter, public :: OQP_xyz_oldcomment = OQP_prefix // "saved geo from previous step" - character(len=*), parameter, public :: all_tags(32) = (/ character(len=80) :: & + character(len=*), parameter, public :: all_tags(36) = (/ character(len=80) :: & OQP_DM_A, OQP_DM_B, OQP_FOCK_A, OQP_FOCK_B, OQP_E_MO_A, OQP_E_MO_B, & OQP_VEC_MO_A, OQP_VEC_MO_B, OQP_Hcore, OQP_SM, OQP_TM, OQP_WAO, & OQP_td_abxc, OQP_td_bvec_mo, OQP_td_mrsf_density, OQP_td_p, OQP_td_t, & OQP_log_filename, OQP_basis_filename, OQP_hbasis_filename, & OQP_xyz_old, OQP_overlap_mo, OQP_overlap_ao, OQP_E_MO_A_old, OQP_E_MO_B_old, & OQP_VEC_MO_A_old, OQP_VEC_MO_B_old, OQP_td_bvec_mo_old, OQP_td_energies_old, & - OQP_nac, OQP_td_states_phase, OQP_td_states_overlap /) + OQP_nac, OQP_td_states_phase, OQP_td_states_overlap, & + OQP_Hqmmm, OQP_mm_potential, OQP_partial_charges,OQP_mm_energy /) interface tagarray_get_data module procedure tagarray_get_data_int64_val, tagarray_get_data_int64_1d, tagarray_get_data_int64_2d, tagarray_get_data_int64_3d module procedure tagarray_get_data_real64_val, tagarray_get_data_real64_1d, tagarray_get_data_real64_2d, tagarray_get_data_real64_3d diff --git a/source/types.F90 b/source/types.F90 index 2f16609..0a9ced3 100644 --- a/source/types.F90 +++ b/source/types.F90 @@ -121,6 +121,7 @@ module types real(c_double) :: resp_constr = 0.01 !< RESP charges constraint logical(c_bool) :: basis_set_issue = .false. !< Basis set issue flag real(c_double) :: conf_print_threshold = 5.0d-02 !< The threshold for configuration printout + logical(c_bool) :: qmmm_flag = .false. !< QM/MM Flag end type control_parameters type, public, bind(c) :: tddft_parameters