From: Adam Sieradzan Date: Tue, 9 Oct 2018 10:43:20 +0000 (+0200) Subject: NEW CORR added to UNRES4 X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=commitdiff_plain;h=817262c16c50ae5848d85d1162e4977c700e1b6f;hp=a515d0fe469dc0bda93e29298c9d981a6e25c603;p=unres4.git NEW CORR added to UNRES4 --- diff --git a/CMakeLists.txt b/CMakeLists.txt index de3836d..de9672a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -114,8 +114,8 @@ set(CMAKE_INSTALL_PREFIX "${CMAKE_SOURCE_DIR}/bin" CACHE PATH "Binary install di # Set force field if (NOT UNRES_FF) - set(UNRES_MD_FF "GAB" CACHE STRING "Choose the force field, options are: GAB E0LL2Y" ) - set_property(CACHE UNRES_MD_FF PROPERTY STRINGS "GAB" "E0LL2Y") + set(UNRES_MD_FF "GAB" CACHE STRING "Choose the force field, options are: GAB E0LL2Y NEWCORR" ) + set_property(CACHE UNRES_MD_FF PROPERTY STRINGS "GAB" "E0LL2Y NEWCORR") endif (NOT UNRES_FF) # Use of MPI library (default ON) diff --git a/PARAM/fourier_opt.parm.OPT_TRP1_FSD_Villin_E0L_QHK_N9L_LX7_BDD_I18 b/PARAM/fourier_opt.parm.OPT_TRP1_FSD_Villin_E0L_QHK_N9L_LX7_BDD_I18 new file mode 100644 index 0000000..4b0e736 --- /dev/null +++ b/PARAM/fourier_opt.parm.OPT_TRP1_FSD_Villin_E0L_QHK_N9L_LX7_BDD_I18 @@ -0,0 +1,219 @@ + -3 # Number of local interaction types + 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 2 1 1 0 1 1 1 1 + 10 9 20 25 +GLY + -0.8621452232 b11(1) + 0.0000000000 b12(1) + 0.7313789248 b11(2) + 0.0000000000 b12(2) + 4.0465105732 b11(3) + 0.0000000000 b12(3) + -3.2270791580 b21(1) + 0.0000000000 b22(1) + -0.7221426339 b21(2) + 0.0000000000 b22(2) + 5.6908471718 b21(3) + 0.0000000000 b22(3) + 10.3828800000 c11(1) + 0.0000000000 c12(1) + 21.0198400000 c11(2) + 0.0000000000 c12(2) + 0.0000000000 c11(3) + 0.0000000000 c12(3) + 0.3056400000 d11(1) + 0.0000000000 d12(1) + -0.2647900000 d11(2) + 0.0000000000 d12(2) + 0.0000000000 d11(3) + 0.0000000000 d12(3) + 1.5149493993 e11(1) + 0.0000000000 e12(1) + 0.0000000000 e21(1) + -6.1446718210 e22(1) + 7.0020624196 e11(2) + 0.0000000000 e12(2) + 0.0000000000 e21(2) + 6.1758090305 e22(2) + -4.2731718047 e0(1) + 0.0000000000 e0(2) + 0.0000000000 e0(3) +ALA + -0.2878386542 b11(1) + 0.5340246908 b12(1) + -0.3425708486 b11(2) + -3.9240673068 b12(2) + 4.2593797716 b11(3) + 1.6082290564 b12(3) + -0.6092499813 b21(1) + 1.1810741056 b22(1) + -2.8679137198 b21(2) + 2.8601711713 b22(2) + 3.6252573043 b21(3) + -13.2722736676 b22(3) + 7.5561100000 c11(1) + 0.5482200000 c12(1) + 12.0045300000 c11(2) + 15.7514800000 c12(2) + 0.0000000000 c11(3) + 0.0000000000 c12(3) + 0.2337300000 d11(1) + 0.0230900000 d12(1) + 0.3581400000 d11(2) + -0.5413800000 d12(2) + 0.0000000000 d11(3) + 0.0000000000 d12(3) + 1.5427601483 e11(1) + -2.0946664525 e12(1) + -3.7118682681 e21(1) + -3.2423931510 e22(1) + 1.1561341780 e11(2) + -2.8637053518 e12(2) + 0.2964859036 e21(2) + 2.0000658209 e22(2) + -1.5289105641 e0(1) + 1.5045800591 e0(2) + 2.1703940418 e0(3) +PRO + -1.5707308981 b11(1) + 0.5008580090 b12(1) + 0.1295058807 b11(2) + 1.1878038803 b12(2) + 4.4917056098 b11(3) + 0.2100906598 b12(3) + 0.7242132444 b21(1) + -3.2698632262 b22(1) + -2.2572254926 b21(2) + 5.2569226538 b22(2) + 5.0417388407 b21(3) + 5.7513919793 b22(3) + 13.9941200000 c11(1) + 5.5338100000 c12(1) + 4.4200400000 c11(2) + -25.4290800000 c12(2) + 0.0000000000 c11(3) + 0.0000000000 c12(3) + 0.6454300000 d11(1) + 0.5700300000 d12(1) + 1.2973700000 d11(2) + 0.7918800000 d12(2) + 0.0000000000 d11(3) + 0.0000000000 d12(3) + 0.5924944193 e11(1) + -3.2523385083 e12(1) + -2.0293467488 e21(1) + 1.3690442042 e22(1) + -0.5531839247 e11(2) + -5.4469248588 e12(2) + -3.5578194175 e21(2) + 0.9317664393 e22(2) + 1.5599919952 e0(1) + 3.1444815253 e0(2) + 2.3755044763 e0(3) +GLY + -0.3365608821 b11(1) + 0.0000000000 b12(1) + -1.7815212118 b11(2) + 0.0000000000 b12(2) + 2.8654373865 b11(3) + 0.0000000000 b12(3) + -1.7125429107 b21(1) + 0.0000000000 b22(1) + -2.2117865414 b21(2) + 0.0000000000 b22(2) + 4.1244646149 b21(3) + 0.0000000000 b22(3) + 0.1769155681 c11(1) + 0.0000000000 c12(1) + 18.3022694952 c11(2) + 0.0000000000 c12(2) + 0.0000000000 c11(3) + 0.0000000000 c12(3) + 0.6479313147 d11(1) + 0.0000000000 d12(1) + 0.3123945472 d11(2) + 0.0000000000 d12(2) + 0.0000000000 d11(3) + 0.0000000000 d12(3) + 3.1692773766 e11(1) + 0.0000000000 e12(1) + 0.0000000000 e21(1) + -4.7442619198 e22(1) + 6.0829185055 e11(2) + 0.0000000000 e12(2) + 0.0000000000 e21(2) + 9.3047137272 e22(2) + -2.4585407292 e0(1) + 0.0000000000 e0(2) + 0.0000000000 e0(3) +ALA + 0.0699903800 b11(1) + 0.6473400960 b12(1) + 0.6307567638 b11(2) + -3.3038484402 b12(2) + 3.6985100887 b11(3) + 0.6456466448 b12(3) + -0.7776661176 b21(1) + 0.1656030542 b22(1) + -1.7179510926 b21(2) + -0.5840593818 b22(2) + 1.6791371601 b21(3) + -0.2407561887 b22(3) + 4.4008957899 c11(1) + -1.3966832985 c12(1) + 5.6452946950 c11(2) + 13.3253230560 c12(2) + 0.0000000000 c11(3) + 0.0000000000 c12(3) + 0.2587893240 d11(1) + 0.1198381815 d12(1) + -0.2151230671 d11(2) + -1.1482082176 d12(2) + 0.0000000000 d11(3) + 0.0000000000 d12(3) + 1.3864046666 e11(1) + 0.0706805962 e12(1) + -1.8747963273 e21(1) + -1.0171241446 e22(1) + 0.6833940808 e11(2) + -3.6214085316 e12(2) + -0.9493445880 e21(2) + 0.6836133612 e22(2) + -0.1136367361 e0(1) + 0.8334023938 e0(2) + 1.2586019646 e0(3) +PRO + -1.2653242255 b11(1) + 4.3252251756 b12(1) + 3.5877938324 b11(2) + 1.9134262096 b12(2) + 3.0640982419 b11(3) + 5.9069671174 b12(3) + 0.2946509372 b21(1) + 0.1132213433 b22(1) + -0.3498120634 b21(2) + -0.2815716594 b22(2) + 1.0143464807 b21(3) + -0.4663010146 b22(3) + 11.5809307622 c11(1) + 7.9390520348 c12(1) + -0.2631081246 c11(2) + -25.3395396435 c12(2) + 0.0000000000 c11(3) + 0.0000000000 c12(3) + 0.4967227582 d11(1) + 0.5273390759 d12(1) + 1.2707336194 d11(2) + 0.5349975067 d12(2) + 0.0000000000 d11(3) + 0.0000000000 d12(3) + 0.3586504760 e11(1) + -2.6779516657 e12(1) + -1.9268413854 e21(1) + 0.8288233035 e22(1) + 0.1311826103 e11(2) + -5.9712643940 e12(2) + -4.9858015847 e21(2) + 0.1769615908 e22(2) + 0.4413749734 e0(1) + 3.0497509063 e0(2) + 2.5125369772 e0(3) diff --git a/PARAM/theta_opt.parm.OPT_TRP1_FSD_Villin_E0L_QHK_N9L_LX7_BDD_I18 b/PARAM/theta_opt.parm.OPT_TRP1_FSD_Villin_E0L_QHK_N9L_LX7_BDD_I18 new file mode 100644 index 0000000..34e6358 --- /dev/null +++ b/PARAM/theta_opt.parm.OPT_TRP1_FSD_Villin_E0L_QHK_N9L_LX7_BDD_I18 @@ -0,0 +1,71 @@ + 3 +12 ************ p + 0 348.2643220864 + 1 653.9849593975 + 2 602.7471290817 + 3 504.6718813684 + 4 409.9349875339 + 5 306.2663644893 + 6 217.1797937290 + 7 140.4770599539 + 8 84.5739144114 + 9 43.7102735331 + 10 20.8302780646 + 11 7.0063070222 + 12 2.4789692685 +12 ************ a + 0 2724.0077280056 + 1 5228.3265790979 + 2 4653.5625905064 + 3 3816.1556806065 + 4 2886.4019534026 + 5 1994.8148094576 + 6 1253.6156070751 + 7 703.6943002523 + 8 348.1657829457 + 9 146.1235296937 + 10 50.3525621590 + 11 12.6691695761 + 12 2.0033525161 +12 ************ G + 0 12981.1654212053 + 1 24905.6646638826 + 2 22008.5042923203 + 3 17862.1792098516 + 4 13281.6538712778 + 5 8989.1340928037 + 6 5496.5892109793 + 7 2993.7693700741 + 8 1428.9884692856 + 9 578.1859975736 + 10 190.5236538627 + 11 46.1013457759 + 12 7.0681037782 +12 ************ A + 0 2724.0077280056 + 1 5228.3265790979 + 2 4653.5625905064 + 3 3816.1556806065 + 4 2886.4019534026 + 5 1994.8148094576 + 6 1253.6156070751 + 7 703.6943002523 + 8 348.1657829457 + 9 146.1235296937 + 10 50.3525621590 + 11 12.6691695761 + 12 2.0033525161 +12 ************ P + 0 348.2643220864 + 1 653.9849593975 + 2 602.7471290817 + 3 504.6718813684 + 4 409.9349875339 + 5 306.2663644893 + 6 217.1797937290 + 7 140.4770599539 + 8 84.5739144114 + 9 43.7102735331 + 10 20.8302780646 + 11 7.0063070222 + 12 2.4789692685 diff --git a/PARAM/torsion_abinitio.parm-2d-all-DL-03-02-2cos b/PARAM/torsion_abinitio.parm-2d-all-DL-03-02-2cos new file mode 100644 index 0000000..a2b3903 --- /dev/null +++ b/PARAM/torsion_abinitio.parm-2d-all-DL-03-02-2cos @@ -0,0 +1,502 @@ +3 +1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 2 1 1 0 1 +************ p-p + 2 3 + 1 0 0 2.46261E+00 -7.99262E-01 + 1 0 1 -4.40162E+00 8.17613E-01 + 1 0 2 -2.91288E+00 3.90340E+00 + 1 1 0 -1.33212E+00 1.48032E-01 + 1 1 1 3.43800E+00 -1.35531E-01 + 1 1 2 1.72784E+00 -3.25189E+00 + 1 2 0 -6.09066E+00 -2.85544E+00 + 1 2 1 1.22093E+01 -6.48650E+00 + 1 2 2 4.23822E-01 -2.79189E+01 + 2 0 0 3.71192E+00 3.06950E+00 + 2 0 1 -1.87483E+00 -3.74230E+00 + 2 0 2 -5.19934E+00 -5.48463E+00 + 2 1 0 5.03153E-01 6.05081E+00 + 2 1 1 -1.83539E+00 -8.20441E+00 + 2 1 2 -9.96238E+00 -9.29723E+00 + 2 2 0 -1.02706E+01 -5.58902E+00 + 2 2 1 1.11689E+01 2.84133E+00 + 2 2 2 2.22836E+01 1.35871E+01 +************ p-a + 2 3 + 1 0 0 1.21730E+00 -1.08644E+00 + 1 0 1 -2.93893E+00 -2.14296E-01 + 1 0 2 -2.91549E+00 9.95310E+00 + 1 1 0 2.40225E+00 1.98200E-01 + 1 1 1 2.38373E-02 -2.75470E-01 + 1 1 2 -3.72367E+00 -6.37079E+00 + 1 2 0 -3.58562E+00 1.46444E+00 + 1 2 1 3.73290E+00 5.73969E+00 + 1 2 2 2.08734E+01 -5.09504E+01 + 2 0 0 -1.81313E-01 1.66637E+00 + 2 0 1 -4.18952E-01 2.97668E+00 + 2 0 2 -5.01059E+00 -1.55791E+00 + 2 1 0 -5.10309E-01 2.43111E+00 + 2 1 1 -1.23136E+00 6.58418E+00 + 2 1 2 -3.36829E+00 -5.44780E+00 + 2 2 0 5.12870E+00 -1.67983E+00 + 2 2 1 7.86879E+00 -2.89554E+00 + 2 2 2 2.40664E+00 -7.12761E+00 +************ p-G + 2 3 + 1 0 0 2.19158E-01 -1.78448E+00 + 1 0 1 -1.06459E-01 2.58308E-01 + 1 0 2 -3.03725E+00 1.18690E+01 + 1 1 0 2.84407E+00 3.91409E-01 + 1 1 1 1.32166E+00 7.09687E-01 + 1 1 2 -4.61750E+00 -1.89442E+00 + 1 2 0 -3.96705E+00 7.04603E+00 + 1 2 1 -4.19641E-01 -7.98555E-01 + 1 2 2 4.37116E+01 -5.54201E+01 + 2 0 0 -1.29673E+00 1.36385E+00 + 2 0 1 -5.29811E-01 2.43939E+00 + 2 0 2 8.83676E-01 9.24388E-01 + 2 1 0 -1.08955E+00 3.29973E+00 + 2 1 1 -6.57979E-01 2.20869E+00 + 2 1 2 2.63344E+00 -5.76594E+00 + 2 2 0 6.51130E+00 1.15217E+00 + 2 2 1 5.26736E+00 -5.65894E+00 + 2 2 2 7.75278E-01 -3.68348E+01 +************ p-A + 2 3 + 1 0 0 -7.29377E-01 -1.32053E+00 + 1 0 1 3.57819E+00 9.54602E-01 + 1 0 2 -1.57774E+00 7.12849E+00 + 1 1 0 4.20580E-01 1.45418E-01 + 1 1 1 -2.22046E+00 4.32348E-01 + 1 1 2 -1.72816E+00 2.87210E+00 + 1 2 0 -1.73035E+00 8.36485E+00 + 1 2 1 -6.77376E+00 -7.30436E+00 + 1 2 2 4.25074E+01 -2.30951E+01 + 2 0 0 -1.96733E+00 1.59714E+00 + 2 0 1 -1.19757E+00 1.31353E+00 + 2 0 2 9.77041E+00 -2.08103E+00 + 2 1 0 -1.31483E+00 1.74297E+00 + 2 1 1 -1.12566E+00 -2.56576E+00 + 2 1 2 5.66654E+00 -4.58523E+00 + 2 2 0 6.90439E+00 -2.15111E+00 + 2 2 1 7.74574E+00 -3.25890E+00 + 2 2 2 -2.03638E+01 -8.70268E-01 +************ p-P + 2 3 + 1 0 0 -1.95197E+00 -2.16515E+00 + 1 0 1 3.25949E+00 5.89160E-01 + 1 0 2 2.54763E+00 -2.32420E+00 + 1 1 0 4.59409E+00 1.04906E+00 + 1 1 1 -5.87236E+00 -1.20104E+00 + 1 1 2 -8.29528E+00 8.53729E+00 + 1 2 0 3.36720E+00 1.05962E+01 + 1 2 1 -8.09452E+00 -6.50731E+00 + 1 2 2 -3.64606E+00 3.42576E+01 + 2 0 0 -4.03809E+00 1.70791E+00 + 2 0 1 1.11770E+00 -2.26112E+00 + 2 0 2 1.53341E+00 -2.38469E+00 + 2 1 0 -7.68847E-01 -2.87836E+00 + 2 1 1 1.55311E+00 4.21802E+00 + 2 1 2 -1.14235E+01 4.75466E+00 + 2 2 0 1.71051E+01 -3.88112E+00 + 2 2 1 1.20249E+00 5.26769E+00 + 2 2 2 -2.17841E+01 3.84505E+00 +************ a-p + 2 3 + 1 0 0 3.83054E+00 2.72616E+00 + 1 0 1 -4.74046E+00 6.20620E-01 + 1 0 2 -6.87506E+00 6.13097E+00 + 1 1 0 1.55984E+00 1.31796E+00 + 1 1 1 -3.56481E+00 -2.19585E+00 + 1 1 2 -2.46352E+00 7.34594E+00 + 1 2 0 -8.22423E+00 -1.77034E+01 + 1 2 1 8.08530E+00 -7.00607E+00 + 1 2 2 2.46508E+01 -1.95884E+01 + 2 0 0 -1.05922E+00 1.64824E+00 + 2 0 1 -4.77034E-01 -2.50084E+00 + 2 0 2 -5.46490E+00 -3.16844E+00 + 2 1 0 7.53080E-01 1.26703E+00 + 2 1 1 -6.56878E-01 -3.33964E+00 + 2 1 2 -3.69700E+00 -8.41962E-01 + 2 2 0 7.48427E+00 -2.21223E+00 + 2 2 1 -4.68664E-02 8.99327E-01 + 2 2 2 5.91888E+00 1.14903E+01 +************ a-a + 2 3 + 1 0 0 1.53561E+00 4.95284E-01 + 1 0 1 4.92733E-01 -1.08071E+00 + 1 0 2 -5.36597E+00 7.63067E+00 + 1 1 0 4.47467E-01 1.12497E+00 + 1 1 1 9.29209E+00 -1.88859E+00 + 1 1 2 -3.56212E+00 -3.87148E+00 + 1 2 0 -4.72525E+00 -5.32980E+00 + 1 2 1 1.95716E+00 7.91566E+00 + 1 2 2 3.20937E+01 -1.41756E+01 + 2 0 0 -1.26374E+00 -3.51061E-01 + 2 0 1 -8.41002E-01 1.46378E+00 + 2 0 2 2.68702E+00 9.08276E-01 + 2 1 0 -2.85896E-01 -1.90090E+00 + 2 1 1 5.53061E-01 -7.25223E-01 + 2 1 2 1.44371E+00 3.69248E+00 + 2 2 0 5.10613E+00 2.52850E+00 + 2 2 1 7.89378E+00 -1.20836E+00 + 2 2 2 -3.03480E+00 -1.02909E+01 +************ a-G + 2 3 + 1 0 0 1.14420E+00 -1.06914E+00 + 1 0 1 3.22255E+00 -1.92253E-01 + 1 0 2 -5.17968E+00 5.37399E+00 + 1 1 0 3.22124E+00 9.73384E-01 + 1 1 1 6.43595E+00 1.06025E+00 + 1 1 2 -7.97103E+00 -8.43302E+00 + 1 2 0 -5.37171E+00 -3.48636E-01 + 1 2 1 -2.15123E+00 3.94834E+00 + 1 2 2 4.11603E+01 -7.89556E-01 + 2 0 0 -1.53627E+00 -7.84139E-01 + 2 0 1 -3.06108E-01 7.21145E-01 + 2 0 2 2.66113E+00 3.49638E-01 + 2 1 0 -1.44009E-02 -2.03085E+00 + 2 1 1 5.92916E-01 -1.87541E+00 + 2 1 2 3.79320E+00 2.91843E+00 + 2 2 0 6.69463E+00 2.81939E+00 + 2 2 1 4.31102E+00 -3.22544E+00 + 2 2 2 1.00906E+01 8.55041E+00 +************ a-A + 2 3 + 1 0 0 -2.38221E-01 -1.66633E+00 + 1 0 1 5.20451E+00 1.04014E+00 + 1 0 2 -2.42033E+00 -4.62892E-01 + 1 1 0 3.40198E+00 2.62420E-01 + 1 1 1 -4.51805E-01 3.49257E+00 + 1 1 2 -9.62634E+00 -8.66714E+00 + 1 2 0 -1.84867E+00 4.93119E+00 + 1 2 1 -6.40319E+00 -1.10575E+00 + 1 2 2 2.90461E+01 1.54289E+01 + 2 0 0 -1.59935E+00 -5.82348E-01 + 2 0 1 -6.62316E-01 -2.63563E+00 + 2 0 2 4.32246E+00 6.32765E-01 + 2 1 0 -1.47054E-01 -1.37702E+00 + 2 1 1 -2.76596E-01 -4.70324E+00 + 2 1 2 7.47587E-01 1.41947E+00 + 2 2 0 7.50788E+00 1.34016E+00 + 2 2 1 5.79825E+00 5.07591E+00 + 2 2 2 -1.19119E+01 9.23976E+00 +************ a-P + 2 3 + 1 0 0 5.63221E-01 -4.22606E+00 + 1 0 1 7.17695E-01 -4.41795E-01 + 1 0 2 2.18620E-02 -4.94417E+00 + 1 1 0 9.99098E+00 -5.93766E-01 + 1 1 1 -9.48298E+00 1.20621E+00 + 1 1 2 -1.55236E+01 -7.39907E+00 + 1 2 0 -6.25326E-01 1.44491E+01 + 1 2 1 -3.77789E+00 7.64695E+00 + 1 2 2 -1.15550E+01 2.50426E+01 + 2 0 0 -8.99742E-01 -1.84769E+00 + 2 0 1 -2.09666E-01 2.23342E+00 + 2 0 2 -4.50273E+00 3.66834E+00 + 2 1 0 7.46585E-01 -1.36173E+00 + 2 1 1 -8.01578E-01 2.12755E-01 + 2 1 2 -2.92761E+00 1.85001E+00 + 2 2 0 5.70168E+00 3.99745E+00 + 2 2 1 -2.29782E+00 -4.90775E+00 + 2 2 2 4.97603E+00 -1.45490E+01 +************ G-p + 2 3 + 1 0 0 4.03909E+00 7.23458E+00 + 1 0 1 -3.46462E+00 7.43451E-01 + 1 0 2 -6.69309E+00 9.73916E+00 + 1 1 0 8.09285E+00 1.11845E+00 + 1 1 1 -7.55494E+00 -1.17643E+00 + 1 1 2 -1.25639E+01 1.17103E+01 + 1 2 0 -7.85853E+00 -3.04240E+01 + 1 2 1 4.36680E+00 -1.09819E+01 + 1 2 2 1.44611E+01 -3.92698E+01 + 2 0 0 -4.26154E-01 1.85592E+00 + 2 0 1 6.16930E-02 -4.09338E+00 + 2 0 2 -4.20745E+00 -3.61304E+00 + 2 1 0 1.09569E+00 -4.92197E+00 + 2 1 1 -5.59509E-01 2.07628E+00 + 2 1 2 9.96455E+00 9.28142E+00 + 2 2 0 -4.62129E-01 -3.98242E+00 + 2 2 1 -7.80550E+00 8.51017E+00 + 2 2 2 -1.42596E+01 1.48946E+01 +************ G-a + 2 3 + 1 0 0 7.49512E-01 8.76692E-01 + 1 0 1 3.89209E+00 -1.62962E+00 + 1 0 2 -4.26454E+00 5.79966E+00 + 1 1 0 1.69928E+00 4.20247E-01 + 1 1 1 3.75319E+00 -2.43701E+00 + 1 1 2 -4.63145E+00 -1.08068E-01 + 1 2 0 -3.47843E+00 -4.56540E+00 + 1 2 1 -6.39033E+00 8.73753E+00 + 1 2 2 3.13593E+01 -2.02261E+01 + 2 0 0 -6.35482E-01 2.82389E-01 + 2 0 1 8.74653E-01 2.21575E+00 + 2 0 2 3.49618E+00 -4.24615E-01 + 2 1 0 1.63160E+00 -8.14674E-01 + 2 1 1 3.42759E+00 -1.99859E+00 + 2 1 2 -2.21850E+00 8.33359E-04 + 2 2 0 4.07601E+00 1.03303E+00 + 2 2 1 2.43192E+00 -2.36801E+00 + 2 2 2 -1.37593E+01 -1.58483E+01 +************ G-G + 2 3 + 1 0 0 1.69533E+00 0.00000E+00 + 1 0 1 5.47750E+00 0.00000E+00 + 1 0 2 -5.29419E+00 0.00000E+00 + 1 1 0 3.40328E+00 0.00000E+00 + 1 1 1 6.95727E+00 0.00000E+00 + 1 1 2 -5.43520E+00 0.00000E+00 + 1 2 0 -5.79872E+00 0.00000E+00 + 1 2 1 -8.83199E+00 0.00000E+00 + 1 2 2 3.65154E+01 0.00000E+00 + 2 0 0 -6.71638E-01 0.00000E+00 + 2 0 1 6.85960E-01 0.00000E+00 + 2 0 2 6.20018E+00 0.00000E+00 + 2 1 0 1.31671E+00 0.00000E+00 + 2 1 1 1.83576E+00 0.00000E+00 + 2 1 2 4.33783E+00 0.00000E+00 + 2 2 0 5.30648E+00 0.00000E+00 + 2 2 1 3.76737E-01 0.00000E+00 + 2 2 2 3.82717E+00 0.00000E+00 +************ G-A + 2 3 + 1 0 0 7.27828E-01 -9.31173E-01 + 1 0 1 3.91718E+00 1.62243E+00 + 1 0 2 -4.26106E+00 -5.46053E+00 + 1 1 0 1.60522E+00 -4.28211E-01 + 1 1 1 3.58698E+00 2.42650E+00 + 1 1 2 -4.65557E+00 1.23172E-01 + 1 2 0 -3.50067E+00 4.81280E+00 + 1 2 1 -6.35805E+00 -8.62036E+00 + 1 2 2 3.18991E+01 1.86756E+01 + 2 0 0 -6.49671E-01 -2.73031E-01 + 2 0 1 8.17379E-01 -2.24272E+00 + 2 0 2 3.44965E+00 4.24126E-01 + 2 1 0 1.63168E+00 8.04187E-01 + 2 1 1 3.33706E+00 2.52574E+00 + 2 1 2 -2.38653E+00 7.73007E-01 + 2 2 0 4.06201E+00 -1.08329E+00 + 2 2 1 3.22974E+00 2.24619E+00 + 2 2 2 -1.26284E+01 1.60739E+01 +************ G-P + 2 3 + 1 0 0 4.00920E+00 -7.22717E+00 + 1 0 1 -3.43061E+00 -7.26845E-01 + 1 0 2 -6.62270E+00 -9.78451E+00 + 1 1 0 8.07187E+00 -1.13027E+00 + 1 1 1 -7.51501E+00 1.17630E+00 + 1 1 2 -1.25097E+01 -1.15688E+01 + 1 2 0 -7.74758E+00 3.03937E+01 + 1 2 1 4.31125E+00 1.06989E+01 + 1 2 2 1.39273E+01 3.97445E+01 + 2 0 0 -4.25031E-01 -1.85987E+00 + 2 0 1 2.23244E-02 4.06427E+00 + 2 0 2 -4.17115E+00 3.68607E+00 + 2 1 0 1.09617E+00 4.91791E+00 + 2 1 1 -5.86927E-01 -2.12222E+00 + 2 1 2 9.84850E+00 -9.17475E+00 + 2 2 0 -5.78311E-01 4.09742E+00 + 2 2 1 -7.21508E+00 -8.42191E+00 + 2 2 2 -1.45640E+01 -1.58664E+01 +************ A-p + 2 3 + 1 0 0 5.78243E-01 4.18378E+00 + 1 0 1 7.10191E-01 3.90179E-01 + 1 0 2 2.05326E-01 5.00525E+00 + 1 1 0 1.04276E+01 6.37466E-01 + 1 1 1 -9.63504E+00 -1.25865E+00 + 1 1 2 -1.61795E+01 7.28220E+00 + 1 2 0 -5.61098E-01 -1.41716E+01 + 1 2 1 -3.87922E+00 -7.25873E+00 + 1 2 2 -1.33148E+01 -2.52719E+01 + 2 0 0 -9.69110E-01 1.86985E+00 + 2 0 1 -1.41997E-01 -2.24943E+00 + 2 0 2 -4.59246E+00 -3.67440E+00 + 2 1 0 7.77849E-01 1.50664E+00 + 2 1 1 -8.19129E-01 -1.17843E-01 + 2 1 2 -2.71538E+00 -2.10930E+00 + 2 2 0 6.25294E+00 -3.96861E+00 + 2 2 1 -3.05071E+00 4.89189E+00 + 2 2 2 5.58213E+00 1.44003E+01 +************ A-a + 2 3 + 1 0 0 -2.13242E-01 1.64040E+00 + 1 0 1 5.20628E+00 -1.02180E+00 + 1 0 2 -2.50370E+00 5.52514E-01 + 1 1 0 3.46681E+00 -2.91525E-01 + 1 1 1 -5.94045E-01 -3.45041E+00 + 1 1 2 -9.80199E+00 8.78027E+00 + 1 2 0 -1.90970E+00 -4.86161E+00 + 1 2 1 -6.40347E+00 8.84078E-01 + 1 2 2 2.94350E+01 -1.56692E+01 + 2 0 0 -1.62437E+00 5.99696E-01 + 2 0 1 -7.01919E-01 2.60602E+00 + 2 0 2 4.43356E+00 -6.53438E-01 + 2 1 0 -1.33117E-01 1.50602E+00 + 2 1 1 -2.88725E-01 4.81905E+00 + 2 1 2 6.55187E-01 -1.52747E+00 + 2 2 0 7.73106E+00 -1.28571E+00 + 2 2 1 6.21523E+00 -5.08066E+00 + 2 2 2 -1.26700E+01 -9.55143E+00 +************ A-G + 2 3 + 1 0 0 1.14055E+00 1.07301E+00 + 1 0 1 3.22985E+00 1.67869E-01 + 1 0 2 -5.05793E+00 -5.47298E+00 + 1 1 0 3.40172E+00 -9.74488E-01 + 1 1 1 6.58088E+00 -1.13736E+00 + 1 1 2 -8.19959E+00 8.33184E+00 + 1 2 0 -5.34883E+00 3.15985E-01 + 1 2 1 -2.25697E+00 -3.73477E+00 + 1 2 2 4.01618E+01 1.55640E+00 + 2 0 0 -1.53964E+00 7.89423E-01 + 2 0 1 -3.57562E-01 -7.11044E-01 + 2 0 2 2.70907E+00 -4.32087E-01 + 2 1 0 1.77831E-02 2.06784E+00 + 2 1 1 5.22943E-01 1.89139E+00 + 2 1 2 3.58707E+00 -2.96251E+00 + 2 2 0 6.72349E+00 -2.90347E+00 + 2 2 1 4.80387E+00 3.22999E+00 + 2 2 2 9.79094E+00 -7.71471E+00 +************ A-A + 2 3 + 1 0 0 1.52932E+00 -4.83000E-01 + 1 0 1 5.26244E-01 1.10003E+00 + 1 0 2 -5.34791E+00 -7.61864E+00 + 1 1 0 4.85415E-01 -1.12317E+00 + 1 1 1 9.57023E+00 1.84944E+00 + 1 1 2 -3.45662E+00 3.96154E+00 + 1 2 0 -4.76155E+00 5.19148E+00 + 1 2 1 1.78206E+00 -8.01928E+00 + 1 2 2 3.19961E+01 1.47297E+01 + 2 0 0 -1.29615E+00 3.54053E-01 + 2 0 1 -9.52931E-01 -1.50334E+00 + 2 0 2 2.68543E+00 -9.34622E-01 + 2 1 0 -2.58346E-01 1.83690E+00 + 2 1 1 4.54076E-01 7.77389E-01 + 2 1 2 1.03020E+00 -3.48662E+00 + 2 2 0 5.28367E+00 -2.61722E+00 + 2 2 1 9.03343E+00 1.28215E+00 + 2 2 2 -2.79174E+00 1.07459E+01 +************ A-P + 2 3 + 1 0 0 3.82682E+00 -2.79811E+00 + 1 0 1 -4.73590E+00 -5.95104E-01 + 1 0 2 -7.07768E+00 -5.98909E+00 + 1 1 0 1.28769E+00 -1.34012E+00 + 1 1 1 -3.58669E+00 2.23107E+00 + 1 1 2 -2.11587E+00 -7.36564E+00 + 1 2 0 -8.27960E+00 1.83390E+01 + 1 2 1 8.14311E+00 6.87076E+00 + 1 2 2 2.65818E+01 1.81265E+01 + 2 0 0 -1.11048E+00 -1.65526E+00 + 2 0 1 -4.95291E-01 2.49076E+00 + 2 0 2 -5.64034E+00 3.20900E+00 + 2 1 0 7.25906E-01 -1.26402E+00 + 2 1 1 -6.21742E-01 3.40131E+00 + 2 1 2 -3.65510E+00 7.79984E-01 + 2 2 0 7.71295E+00 2.29684E+00 + 2 2 1 1.73800E-01 -8.07209E-01 + 2 2 2 7.65091E+00 -1.18487E+01 +************ P-p + 2 3 + 1 0 0 -2.30263E+00 8.20545E-01 + 1 0 1 3.41741E+00 -2.61345E-01 + 1 0 2 3.14939E+00 1.94187E+00 + 1 1 0 4.33856E+00 -1.09306E+00 + 1 1 1 -5.60761E+00 1.29668E+00 + 1 1 2 -7.89092E+00 -9.66713E+00 + 1 2 0 3.77073E+00 -9.68414E+00 + 1 2 1 -8.28126E+00 6.58108E+00 + 1 2 2 -2.74988E+00 -3.55473E+01 + 2 0 0 -5.16150E+00 -1.08766E+00 + 2 0 1 1.40792E+00 1.78462E+00 + 2 0 2 -3.00686E-01 8.32415E-01 + 2 1 0 -6.27266E-01 3.12325E+00 + 2 1 1 1.70665E+00 -4.94984E+00 + 2 1 2 -1.34932E+01 -5.16953E+00 + 2 2 0 1.86347E+01 3.14407E+00 + 2 2 1 1.83665E+00 -5.02774E+00 + 2 2 2 -2.20926E+01 -2.82096E+00 +************ P-a + 2 3 + 1 0 0 -8.64645E-01 1.12361E+00 + 1 0 1 2.82003E+00 -9.17028E-01 + 1 0 2 -1.16970E+00 -7.58888E+00 + 1 1 0 5.17823E-01 -8.45405E-02 + 1 1 1 -2.10936E+00 -6.24806E-01 + 1 1 2 -1.91327E+00 -3.01746E+00 + 1 2 0 -1.49705E+00 -8.01488E+00 + 1 2 1 -5.34343E+00 7.32081E+00 + 1 2 2 4.19203E+01 2.38729E+01 + 2 0 0 -2.80786E+00 -1.86496E+00 + 2 0 1 -1.73439E+00 -8.33533E-01 + 2 0 2 1.23011E+01 2.78226E+00 + 2 1 0 -1.23430E+00 -1.77876E+00 + 2 1 1 -1.05706E+00 2.39243E+00 + 2 1 2 5.38227E+00 4.59709E+00 + 2 2 0 8.77093E+00 2.84851E+00 + 2 2 1 9.28050E+00 2.34959E+00 + 2 2 2 -2.54933E+01 -1.42375E+00 +************ P-G + 2 3 + 1 0 0 -1.41794E-01 2.24292E+00 + 1 0 1 -8.68328E-01 -2.64070E-02 + 1 0 2 -2.39054E+00 -1.29686E+01 + 1 1 0 2.91533E+00 -3.66684E-01 + 1 1 1 1.69654E+00 -5.93579E-01 + 1 1 2 -4.74338E+00 2.02009E+00 + 1 2 0 -3.30414E+00 -7.92046E+00 + 1 2 1 1.06223E+00 2.31735E-01 + 1 2 2 4.25030E+01 5.75587E+01 + 2 0 0 -2.00147E+00 -1.94996E+00 + 2 0 1 -9.24292E-01 -2.47711E+00 + 2 0 2 2.33877E+00 7.45714E-02 + 2 1 0 -1.06711E+00 -3.22569E+00 + 2 1 1 -8.86338E-01 -2.48453E+00 + 2 1 2 2.32465E+00 5.47203E+00 + 2 2 0 7.93555E+00 5.60872E-02 + 2 2 1 6.28567E+00 5.89176E+00 + 2 2 2 -1.29164E+00 3.54005E+01 +************ P-A + 2 3 + 1 0 0 1.13349E+00 1.27674E+00 + 1 0 1 -3.54991E+00 -5.47760E-02 + 1 0 2 -2.61499E+00 -9.21893E+00 + 1 1 0 2.37580E+00 -2.40150E-01 + 1 1 1 2.56130E-02 5.12615E-01 + 1 1 2 -3.80586E+00 6.24784E+00 + 1 2 0 -3.45701E+00 -1.65685E+00 + 1 2 1 4.93680E+00 -5.17248E+00 + 1 2 2 2.06881E+01 4.92350E+01 + 2 0 0 -7.89653E-01 -2.16183E+00 + 2 0 1 -8.91731E-01 -4.38378E+00 + 2 0 2 -4.02851E+00 2.49568E+00 + 2 1 0 -5.31740E-01 -2.20017E+00 + 2 1 1 -1.33923E+00 -6.29649E+00 + 2 1 2 -3.15716E+00 4.64913E+00 + 2 2 0 6.34569E+00 2.77557E+00 + 2 2 1 9.02181E+00 6.27248E+00 + 2 2 2 1.59213E-01 5.77651E+00 +************ P-P + 2 3 + 1 0 0 1.70628E+00 1.63023E+00 + 1 0 1 -3.73537E+00 -8.13591E-01 + 1 0 2 -1.48057E+00 -3.38357E+00 + 1 1 0 -1.02176E+00 7.03468E-02 + 1 1 1 3.11735E+00 2.56109E-01 + 1 1 2 1.25360E+00 2.86580E+00 + 1 2 0 -4.69677E+00 2.04322E+00 + 1 2 1 1.09575E+01 6.61697E+00 + 1 2 2 -3.64248E+00 2.72021E+01 + 2 0 0 3.60083E+00 -4.31812E+00 + 2 0 1 -2.28071E+00 5.23398E+00 + 2 0 2 -7.13839E+00 7.82677E+00 + 2 1 0 1.43190E-01 -5.70656E+00 + 2 1 1 -1.72594E+00 8.00730E+00 + 2 1 2 -9.81855E+00 8.53652E+00 + 2 2 0 -1.09903E+01 8.44503E+00 + 2 2 1 1.15417E+01 -5.59075E+00 + 2 2 2 2.70959E+01 -2.08570E+01 diff --git a/source/unres/CMakeLists.txt b/source/unres/CMakeLists.txt index 98ec447..ee2fb79 100644 --- a/source/unres/CMakeLists.txt +++ b/source/unres/CMakeLists.txt @@ -111,6 +111,10 @@ if(UNRES_MD_FF STREQUAL "GAB" ) elseif(UNRES_MD_FF STREQUAL "E0LL2Y") # set preprocesor flags set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0" ) +elseif(UNRES_MD_FF STREQUAL "NEWCORR") + # set preprocesor flags + set(CPPFLAGS "PROCOR -DUNRES -DISNAN -DSPLITELE -DLANG0 -DNEWCORR -DCORRCD" ) + #========================================= # Settings for 4P force field diff --git a/source/unres/data/control_data.F90 b/source/unres/data/control_data.F90 index e084359..e4fb4e9 100644 --- a/source/unres/data/control_data.F90 +++ b/source/unres/data/control_data.F90 @@ -30,11 +30,12 @@ ! common /cntrl/ integer :: modecalc,iscode,indpdb,indback,indphi,iranconf,& icheckgrad,iprint,i2ndstr,mucadyn,constr_dist,symetr,shield_mode,& - tubemode,genconstr,afmlog,selfguide,scelemode + tubemode,genconstr,afmlog,selfguide,scelemode,tor_mode logical :: minim,refstr,pdbref,overlapsc,& energy_dec,sideadd,lsecondary,read_cart,unres_pdb,& vdisulf,searchsc,lmuca,dccart,extconf,out1file,& - gnorm_check,gradout,split_ene,with_theta_constr,protein,ions,nucleic + gnorm_check,gradout,split_ene,with_theta_constr,protein,ions,nucleic,& + raw_psipred #ifdef CLUSTER integer :: iopt,nend,nstart,outpdb,outmol2 !cluster logical :: punch_dist,print_dist,lside,lprint_cart,lprint_int,& diff --git a/source/unres/data/energy_data.F90 b/source/unres/data/energy_data.F90 index 5f340cf..6459849 100644 --- a/source/unres/data/energy_data.F90 +++ b/source/unres/data/energy_data.F90 @@ -18,6 +18,11 @@ !----------------------------------------------------------------------------- ! Max. number of SC contacts integer :: maxcont +! Maximum number of valence and torsional in rigorous approach + integer,parameter :: maxtor_kcc=6 + integer,parameter :: maxval_kcc=6 + integer,parameter :: maxang_kcc=36 + !----------------------------------------------------------------------------- ! commom.contacts ! common /contacts/ @@ -68,7 +73,7 @@ wturn6,wvdwpp,wliptran,wshield,lipscale,wtube, & wbond_nucl,wang_nucl,wcorr_nucl,wcorr3_nucl,welpp,wtor_nucl,& wtor_d_nucl,welsb,wsbloc,wvdwsb,welpsb,wvdwpp_nucl,wvdwpsb,wcatprot,& - wcatcat,wscbase,wpepbase,wscpho,wpeppho + wcatcat,wscbase,wpepbase,wscpho,wpeppho,wdihc #ifdef CLUSTER real(kind=8) :: scalscp #endif @@ -295,7 +300,14 @@ real(kind=8),dimension(:,:,:),allocatable :: vlor2,vlor3 !(maxlor,maxtor,maxtor) integer,dimension(:),allocatable :: itortyp !(-ntyp1:ntyp1) integer,dimension(:,:,:),allocatable :: nterm,nlor !(-maxtor:maxtor,-maxtor:maxtor,2) +! ---- for rigorous approach integer :: ntortyp,nterm_old +! integer nloctyp + integer,dimension(:,:),allocatable :: nterm_kcc_Tb,nterm_kcc + integer,dimension(:),allocatable :: iloctyp,itype2loc + real(kind=8),dimension(:,:,:,:,:),allocatable :: v1_kcc,v2_kcc + real(kind=8),dimension(:,:),allocatable :: v1bend_chyb + integer,dimension(:),allocatable :: nbend_kcc_Tb !------torsion nucleic real(kind=8),dimension(:,:),allocatable :: v0_nucl !(-maxtor:maxtor,-maxtor:maxtor,2) real(kind=8),dimension(:,:,:),allocatable :: v1_nucl,v2_nucl !(maxterm,-maxtor:maxtor,-maxtor:maxtor) @@ -317,9 +329,13 @@ ! surfacecommon ! common/fourier/ real(kind=8),dimension(:,:),allocatable :: b1,b2,& - b1tilde !(2,-maxtor:maxtor), + b1tilde,b2tilde,gtb1,gtb2!(2,-maxtor:maxtor), real(kind=8),dimension(:,:,:),allocatable :: cc,dd,ee,& - ctilde,dtilde !(2,2,-maxtor:maxtor) + ctilde,dtilde,bnew1,bnew2,ccnew,ddnew,bnew1tor,& + bnew2tor,ccnewtor,ddnewtor,ccold,ddold,eeold,& + gtCC,gtDD,gtEE,gtEUg + real(kind=8),dimension(:,:,:,:),allocatable :: eenew,eenewtor + real(kind=8),dimension(:,:),allocatable :: e0new,e0newtor integer :: nloctyp ! common/fourier/ z wham real(kind=8),dimension(:,:),allocatable :: b !(13,0:maxtor) @@ -379,4 +395,9 @@ real(kind=8) ,dimension(2) :: wqdip_peppho real(kind=8) :: eps_peppho,sigma_peppho,sigmap1_peppho,sigmap2_peppho, & aa_peppho,bb_peppho +!------------- for psi prec constraints + real(kind=8),dimension(:,:), allocatable :: vpsipred,sdihed + + + end module energy_data diff --git a/source/unres/energy.F90 b/source/unres/energy.F90 index 960b6fe..4e043fe 100644 --- a/source/unres/energy.F90 +++ b/source/unres/energy.F90 @@ -74,7 +74,7 @@ ! amino-acid residue. ! common /precomp1/ real(kind=8),dimension(:,:),allocatable :: mu,muder,Ub2,Ub2der,& - Ctobr,Ctobrder,Dtobr2,Dtobr2der !(2,maxres) + Ctobr,Ctobrder,Dtobr2,Dtobr2der,gUb2 !(2,maxres) real(kind=8),dimension(:,:,:),allocatable :: EUg,EUgder,CUg,& CUgder,DUg,Dugder,DtUg2,DtUg2der !(2,2,maxres) ! This common block contains vectors and matrices dependent on two @@ -87,6 +87,7 @@ real(kind=8),dimension(:,:,:,:),allocatable :: Ug2DtEUgder,& DtUg2EUgder !(2,2,2,maxres) ! common /rotat_old/ + real(kind=8),dimension(4) :: gmuij,gmuij1,gmuij2,gmuji1,gmuji2 real(kind=8),dimension(:),allocatable :: costab,sintab,& costab2,sintab2 !(maxres) ! This common block contains dipole-interaction matrices and their @@ -609,7 +610,7 @@ .or. wcorr4.gt.0.0d0 .or. wcorr5.gt.0.d0 & .or. wcorr6.gt.0.0d0 .or. wturn6.gt.0.0d0 ) then #endif - write(iout,*),"just befor eelec call" +! write(iout,*),"just befor eelec call" call eelec(ees,evdw1,eel_loc,eello_turn3,eello_turn4) ! write (iout,*) "ELEC calc" else @@ -662,12 +663,26 @@ ! Calculate the virtual-bond-angle energy. ! write(iout,*) "in etotal afer edis",ipot - if (wang.gt.0.0d0) then - call ebend(ebe,ethetacnstr) +! if (wang.gt.0.0d0) then +! call ebend(ebe,ethetacnstr) +! else +! ebe=0 +! ethetacnstr=0 +! endif + if (wang.gt.0d0) then + if (tor_mode.eq.0) then + call ebend(ebe) + else +!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the +!C energy function + call ebend_kcc(ebe) + endif else - ebe=0 - ethetacnstr=0 + ebe=0.0d0 endif + ethetacnstr=0.0d0 + if (with_theta_constr) call etheta_constr(ethetacnstr) + ! write(iout,*) "in etotal afer ebe",ipot ! print *,"Processor",myrank," computed UB" @@ -681,12 +696,27 @@ ! Calculate the virtual-bond torsional energy. ! !d print *,'nterm=',nterm - if (wtor.gt.0) then - call etor(etors,edihcnstr) +! if (wtor.gt.0) then +! call etor(etors,edihcnstr) +! else +! etors=0 +! edihcnstr=0 +! endif + if (wtor.gt.0.0d0) then + if (tor_mode.eq.0) then + call etor(etors) + else +!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the +!C energy function + call etor_kcc(etors) + endif else - etors=0 - edihcnstr=0 + etors=0.0d0 endif + edihcnstr=0.0d0 + if (ndih_constr.gt.0) call etor_constr(edihcnstr) +!c print *,"Processor",myrank," computed Utor" + ! print *,"Processor",myrank," computed Utor" ! @@ -2689,7 +2719,8 @@ ! include 'COMMON.FFIELD' real(kind=8) :: auxvec(2),auxmat(2,2) integer :: i,iti1,iti,k,l - real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2 + real(kind=8) :: sin1,cos1,sin2,cos2,dwacos2,dwasin2,cost1,sint1,& + sint1sq,sint1cub,sint1cost1,b1k,b2k,aux ! print *,"in set matrices" ! ! Compute the virtual-bond-torsional-angle dependent quantities needed @@ -2701,6 +2732,149 @@ #else do i=3,nres+1 #endif + if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itype2loc(itype(i-2,1)) + else + iti=nloctyp + endif +!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itype2loc(itype(i-1,1)) + else + iti1=nloctyp + endif +! print *,i,itype(i-2,1),iti +#ifdef NEWCORR + cost1=dcos(theta(i-1)) + sint1=dsin(theta(i-1)) + sint1sq=sint1*sint1 + sint1cub=sint1sq*sint1 + sint1cost1=2*sint1*cost1 +! print *,"cost1",cost1,theta(i-1) +!c write (iout,*) "bnew1",i,iti +!c write (iout,*) (bnew1(k,1,iti),k=1,3) +!c write (iout,*) (bnew1(k,2,iti),k=1,3) +!c write (iout,*) "bnew2",i,iti +!c write (iout,*) (bnew2(k,1,iti),k=1,3) +!c write (iout,*) (bnew2(k,2,iti),k=1,3) + k=1 +! print *,bnew1(1,k,iti),"bnew1" + do k=1,2 + b1k=bnew1(1,k,iti)+(bnew1(2,k,iti)+bnew1(3,k,iti)*cost1)*cost1 +! print *,b1k +! write(*,*) shape(b1) +! if(.not.allocated(b1)) print *, "WTF?" + b1(k,i-2)=sint1*b1k +! +! print *,b1(k,i-2) + + gtb1(k,i-2)=cost1*b1k-sint1sq*& + (bnew1(2,k,iti)+2*bnew1(3,k,iti)*cost1) +! print *,gtb1(k,i-2) + + b2k=bnew2(1,k,iti)+(bnew2(2,k,iti)+bnew2(3,k,iti)*cost1)*cost1 + b2(k,i-2)=sint1*b2k +! print *,b2(k,i-2) + + gtb2(k,i-2)=cost1*b2k-sint1sq*& + (bnew2(2,k,iti)+2*bnew2(3,k,iti)*cost1) +! print *,gtb2(k,i-2) + + enddo +! print *,b1k,b2k + do k=1,2 + aux=ccnew(1,k,iti)+(ccnew(2,k,iti)+ccnew(3,k,iti)*cost1)*cost1 + cc(1,k,i-2)=sint1sq*aux + gtcc(1,k,i-2)=sint1cost1*aux-sint1cub*& + (ccnew(2,k,iti)+2*ccnew(3,k,iti)*cost1) + aux=ddnew(1,k,iti)+(ddnew(2,k,iti)+ddnew(3,k,iti)*cost1)*cost1 + dd(1,k,i-2)=sint1sq*aux + gtdd(1,k,i-2)=sint1cost1*aux-sint1cub*& + (ddnew(2,k,iti)+2*ddnew(3,k,iti)*cost1) + enddo +! print *,"after cc" + cc(2,1,i-2)=cc(1,2,i-2) + cc(2,2,i-2)=-cc(1,1,i-2) + gtcc(2,1,i-2)=gtcc(1,2,i-2) + gtcc(2,2,i-2)=-gtcc(1,1,i-2) + dd(2,1,i-2)=dd(1,2,i-2) + dd(2,2,i-2)=-dd(1,1,i-2) + gtdd(2,1,i-2)=gtdd(1,2,i-2) + gtdd(2,2,i-2)=-gtdd(1,1,i-2) +! print *,"after dd" + + do k=1,2 + do l=1,2 + aux=eenew(1,l,k,iti)+eenew(2,l,k,iti)*cost1 + EE(l,k,i-2)=sint1sq*aux + gtEE(l,k,i-2)=sint1cost1*aux-sint1cub*eenew(2,l,k,iti) + enddo + enddo + EE(1,1,i-2)=EE(1,1,i-2)+e0new(1,iti)*cost1 + EE(1,2,i-2)=EE(1,2,i-2)+e0new(2,iti)+e0new(3,iti)*cost1 + EE(2,1,i-2)=EE(2,1,i-2)+e0new(2,iti)*cost1+e0new(3,iti) + EE(2,2,i-2)=EE(2,2,i-2)-e0new(1,iti) + gtEE(1,1,i-2)=gtEE(1,1,i-2)-e0new(1,iti)*sint1 + gtEE(1,2,i-2)=gtEE(1,2,i-2)-e0new(3,iti)*sint1 + gtEE(2,1,i-2)=gtEE(2,1,i-2)-e0new(2,iti)*sint1 +! print *,"after ee" + +!c b1tilde(1,i-2)=b1(1,i-2) +!c b1tilde(2,i-2)=-b1(2,i-2) +!c b2tilde(1,i-2)=b2(1,i-2) +!c b2tilde(2,i-2)=-b2(2,i-2) +#ifdef DEBUG + write (iout,*) 'i=',i-2,gtb1(2,i-2),gtb1(1,i-2) + write(iout,*) 'b1=',(b1(k,i-2),k=1,2) + write(iout,*) 'b2=',(b2(k,i-2),k=1,2) + write (iout,*) 'theta=', theta(i-1) +#endif +#else + if (i.gt. nnt+2 .and. i.lt.nct+2) then + iti = itype2loc(itype(i-2,1)) + else + iti=nloctyp + endif +!c write (iout,*) "i",i-1," itype",itype(i-2)," iti",iti +!c if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then + if (i.gt. nnt+1 .and. i.lt.nct+1) then + iti1 = itype2loc(itype(i-1,1)) + else + iti1=nloctyp + endif + b1(1,i-2)=b(3,iti) + b1(2,i-2)=b(5,iti) + b2(1,i-2)=b(2,iti) + b2(2,i-2)=b(4,iti) + do k=1,2 + do l=1,2 + CC(k,l,i-2)=ccold(k,l,iti) + DD(k,l,i-2)=ddold(k,l,iti) + EE(k,l,i-2)=eeold(k,l,iti) + enddo + enddo +#endif + b1tilde(1,i-2)= b1(1,i-2) + b1tilde(2,i-2)=-b1(2,i-2) + b2tilde(1,i-2)= b2(1,i-2) + b2tilde(2,i-2)=-b2(2,i-2) +!c + Ctilde(1,1,i-2)= CC(1,1,i-2) + Ctilde(1,2,i-2)= CC(1,2,i-2) + Ctilde(2,1,i-2)=-CC(2,1,i-2) + Ctilde(2,2,i-2)=-CC(2,2,i-2) +!c + Dtilde(1,1,i-2)= DD(1,1,i-2) + Dtilde(1,2,i-2)= DD(1,2,i-2) + Dtilde(2,1,i-2)=-DD(2,1,i-2) + Dtilde(2,2,i-2)=-DD(2,2,i-2) + enddo +#ifdef PARMAT + do i=ivec_start+2,ivec_end+2 +#else + do i=3,nres+1 +#endif + ! print *,i,"i" if (i .lt. nres+1) then sin1=dsin(phi(i)) @@ -2773,37 +2947,43 @@ if (itype(i-2,1).eq.0) then iti=ntortyp+1 else - iti = itortyp(itype(i-2,1)) + iti = itype2loc(itype(i-2,1)) endif else - iti=ntortyp+1 + iti=nloctyp endif ! if (i.gt. iatel_s+1 .and. i.lt.iatel_e+4) then if (i.gt. nnt+1 .and. i.lt.nct+1) then if (itype(i-1,1).eq.0) then - iti1=ntortyp+1 + iti1=nloctyp else - iti1 = itortyp(itype(i-1,1)) + iti1 = itype2loc(itype(i-1,1)) endif else - iti1=ntortyp+1 + iti1=nloctyp endif ! print *,iti,i,"iti",iti1,itype(i-1,1),itype(i-2,1) !d write (iout,*) '*******i',i,' iti1',iti -!d write (iout,*) 'b1',b1(:,iti) -!d write (iout,*) 'b2',b2(:,iti) +! write (iout,*) 'b1',b1(:,iti) +! write (iout,*) 'b2',b2(:,i-2) !d write (iout,*) 'Ug',Ug(:,:,i-2) ! if (i .gt. iatel_s+2) then if (i .gt. nnt+2) then - call matvec2(Ug(1,1,i-2),b2(1,iti),Ub2(1,i-2)) - call matmat2(EE(1,1,iti),Ug(1,1,i-2),EUg(1,1,i-2)) + call matvec2(Ug(1,1,i-2),b2(1,i-2),Ub2(1,i-2)) +#ifdef NEWCORR + call matvec2(Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2)) +!c write (iout,*) Ug(1,1,i-2),gtb2(1,i-2),gUb2(1,i-2),"chuj" +#endif + + call matmat2(EE(1,1,i-2),Ug(1,1,i-2),EUg(1,1,i-2)) + call matmat2(gtEE(1,1,i-2),Ug(1,1,i-2),gtEUg(1,1,i-2)) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or. wcorr6.gt.0.0d0) & then - call matmat2(CC(1,1,iti),Ug(1,1,i-2),CUg(1,1,i-2)) - call matmat2(DD(1,1,iti),Ug(1,1,i-2),DUg(1,1,i-2)) - call matmat2(Dtilde(1,1,iti),Ug2(1,1,i-2),DtUg2(1,1,i-2)) - call matvec2(Ctilde(1,1,iti1),obrot(1,i-2),Ctobr(1,i-2)) - call matvec2(Dtilde(1,1,iti),obrot2(1,i-2),Dtobr2(1,i-2)) + call matmat2(CC(1,1,i-2),Ug(1,1,i-2),CUg(1,1,i-2)) + call matmat2(DD(1,1,i-2),Ug(1,1,i-2),DUg(1,1,i-2)) + call matmat2(Dtilde(1,1,i-2),Ug2(1,1,i-2),DtUg2(1,1,i-2)) + call matvec2(Ctilde(1,1,i-1),obrot(1,i-2),Ctobr(1,i-2)) + call matvec2(Dtilde(1,1,i-2),obrot2(1,i-2),Dtobr2(1,i-2)) endif else do k=1,2 @@ -2818,8 +2998,8 @@ enddo enddo endif - call matvec2(Ugder(1,1,i-2),b2(1,iti),Ub2der(1,i-2)) - call matmat2(EE(1,1,iti),Ugder(1,1,i-2),EUgder(1,1,i-2)) + call matvec2(Ugder(1,1,i-2),b2(1,i-2),Ub2der(1,i-2)) + call matmat2(EE(1,1,i-2),Ugder(1,1,i-2),EUgder(1,1,i-2)) do k=1,2 muder(k,i-2)=Ub2der(k,i-2) enddo @@ -2828,34 +3008,34 @@ if (itype(i-1,1).eq.0) then iti1=ntortyp+1 elseif (itype(i-1,1).le.ntyp) then - iti1 = itortyp(itype(i-1,1)) + iti1 = itype2loc(itype(i-1,1)) else - iti1=ntortyp+1 + iti1=nloctyp endif else - iti1=ntortyp+1 + iti1=nloctyp endif do k=1,2 - mu(k,i-2)=Ub2(k,i-2)+b1(k,iti1) + mu(k,i-2)=Ub2(k,i-2)+b1(k,i-1) enddo -! if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2) -! if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,iti1) -! if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2) + if (energy_dec) write (iout,*) 'Ub2 ',i,Ub2(:,i-2) + if (energy_dec) write (iout,*) 'b1 ',iti1,b1(:,i-1) + if (energy_dec) write (iout,*) 'mu ',i,iti1,mu(:,i-2) !d write (iout,*) 'mu1',mu1(:,i-2) !d write (iout,*) 'mu2',mu2(:,i-2) if (wcorr4.gt.0.0d0 .or. wcorr5.gt.0.0d0 .or.wcorr6.gt.0.0d0) & then - call matmat2(CC(1,1,iti1),Ugder(1,1,i-2),CUgder(1,1,i-2)) - call matmat2(DD(1,1,iti),Ugder(1,1,i-2),DUgder(1,1,i-2)) - call matmat2(Dtilde(1,1,iti),Ug2der(1,1,i-2),DtUg2der(1,1,i-2)) - call matvec2(Ctilde(1,1,iti1),obrot_der(1,i-2),Ctobrder(1,i-2)) - call matvec2(Dtilde(1,1,iti),obrot2_der(1,i-2),Dtobr2der(1,i-2)) + call matmat2(CC(1,1,i-2),Ugder(1,1,i-2),CUgder(1,1,i-2)) + call matmat2(DD(1,1,i-2),Ugder(1,1,i-2),DUgder(1,1,i-2)) + call matmat2(Dtilde(1,1,i-2),Ug2der(1,1,i-2),DtUg2der(1,1,i-2)) + call matvec2(Ctilde(1,1,i-2),obrot_der(1,i-2),Ctobrder(1,i-2)) + call matvec2(Dtilde(1,1,i-2),obrot2_der(1,i-2),Dtobr2der(1,i-2)) ! Vectors and matrices dependent on a single virtual-bond dihedral. - call matvec2(DD(1,1,iti),b1tilde(1,iti1),auxvec(1)) + call matvec2(DD(1,1,i-2),b1tilde(1,iti1),auxvec(1)) call matvec2(Ug2(1,1,i-2),auxvec(1),Ug2Db1t(1,i-2)) call matvec2(Ug2der(1,1,i-2),auxvec(1),Ug2Db1tder(1,i-2)) - call matvec2(CC(1,1,iti1),Ub2(1,i-2),CUgb2(1,i-2)) - call matvec2(CC(1,1,iti1),Ub2der(1,i-2),CUgb2der(1,i-2)) + call matvec2(CC(1,1,i-2),Ub2(1,i-2),CUgb2(1,i-2)) + call matvec2(CC(1,1,i-2),Ub2der(1,i-2),CUgb2der(1,i-2)) call matmat2(EUg(1,1,i-2),CC(1,1,iti1),EUgC(1,1,i-2)) call matmat2(EUgder(1,1,i-2),CC(1,1,iti1),EUgCder(1,1,i-2)) call matmat2(EUg(1,1,i-2),DD(1,1,iti1),EUgD(1,1,i-2)) @@ -3476,6 +3656,7 @@ real(kind=8),dimension(2,2) :: acipa !el,a_temp !el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1 real(kind=8),dimension(4) :: muij + real(kind=8) :: geel_loc_ij,geel_loc_ji real(kind=8) :: xj_safe,yj_safe,zj_safe,xj_temp,yj_temp,zj_temp,& dist_temp, dist_init,rlocshield,fracinbuf integer xshift,yshift,zshift,ilist,iresshield @@ -3917,6 +4098,15 @@ do l=1,2 kkk=kkk+1 muij(kkk)=mu(k,i)*mu(l,j) +#ifdef NEWCORR + gmuij1(kkk)=gtb1(k,i+1)*mu(l,j) +!c write(iout,*) 'k=',k,i,gtb1(k,i+1),gtb1(k,i+1)*mu(l,j) + gmuij2(kkk)=gUb2(k,i)*mu(l,j) + gmuji1(kkk)=mu(k,i)*gtb1(l,j+1) +!c write(iout,*) 'l=',l,j,gtb1(l,j+1),gtb1(l,j+1)*mu(k,i) + gmuji2(kkk)=mu(k,i)*gUb2(l,j) +#endif + enddo enddo !d write (iout,*) 'EELEC: i',i,' j',j @@ -4143,6 +4333,54 @@ enddo endif +#ifdef NEWCORR + geel_loc_ij=(a22*gmuij1(1)& + +a23*gmuij1(2)& + +a32*gmuij1(3)& + +a33*gmuij1(4))& + *fac_shield(i)*fac_shield(j) +!c write(iout,*) "derivative over thatai" +!c write(iout,*) a22*gmuij1(1), a23*gmuij1(2) ,a32*gmuij1(3), +!c & a33*gmuij1(4) + gloc(nphi+i,icg)=gloc(nphi+i,icg)+& + geel_loc_ij*wel_loc +!c write(iout,*) "derivative over thatai-1" +!c write(iout,*) a22*gmuij2(1), a23*gmuij2(2) ,a32*gmuij2(3), +!c & a33*gmuij2(4) + geel_loc_ij=& + a22*gmuij2(1)& + +a23*gmuij2(2)& + +a32*gmuij2(3)& + +a33*gmuij2(4) + gloc(nphi+i-1,icg)=gloc(nphi+i-1,icg)+& + geel_loc_ij*wel_loc& + *fac_shield(i)*fac_shield(j) + +!c Derivative over j residue + geel_loc_ji=a22*gmuji1(1)& + +a23*gmuji1(2)& + +a32*gmuji1(3)& + +a33*gmuji1(4) +!c write(iout,*) "derivative over thataj" +!c write(iout,*) a22*gmuji1(1), a23*gmuji1(2) ,a32*gmuji1(3), +!c & a33*gmuji1(4) + + gloc(nphi+j,icg)=gloc(nphi+j,icg)+& + geel_loc_ji*wel_loc& + *fac_shield(i)*fac_shield(j) + + geel_loc_ji=& + +a22*gmuji2(1)& + +a23*gmuji2(2)& + +a32*gmuji2(3)& + +a33*gmuji2(4) +!c write(iout,*) "derivative over thataj-1" +!c write(iout,*) a22*gmuji2(1), a23*gmuji2(2) ,a32*gmuji2(3), +!c & a33*gmuji2(4) + gloc(nphi+j-1,icg)=gloc(nphi+j-1,icg)+& + geel_loc_ji*wel_loc& + *fac_shield(i)*fac_shield(j) +#endif ! write (iout,*) 'i',i,' j',j,' eel_loc_ij',eel_loc_ij ! eel_loc_ij=0.0 @@ -4501,7 +4739,9 @@ ! include 'COMMON.CONTROL' real(kind=8),dimension(3) :: ggg real(kind=8),dimension(2,2) :: auxmat,auxmat1,auxmat2,pizda,& - e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2 + e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2,gpizda1,& + gpizda2,auxgmat1,auxgmatt1,auxgmat2,auxgmatt2 + real(kind=8),dimension(2) :: auxvec,auxvec1 !el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1 real(kind=8),dimension(2,2) :: auxmat3 !el, a_temp @@ -4561,8 +4801,15 @@ !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !d call checkint_turn3(i,a_temp,eello_turn3_num) call matmat2(EUg(1,1,i+1),EUg(1,1,i+2),auxmat(1,1)) + call matmat2(gtEUg(1,1,i+1),EUg(1,1,i+2),auxgmat1(1,1)) + call matmat2(EUg(1,1,i+1),gtEUg(1,1,i+2),auxgmat2(1,1)) call transpose2(auxmat(1,1),auxmat1(1,1)) + call transpose2(auxgmat1(1,1),auxgmatt1(1,1)) + call transpose2(auxgmat2(1,1),auxgmatt2(1,1)) call matmat2(a_temp(1,1),auxmat1(1,1),pizda(1,1)) + call matmat2(a_temp(1,1),auxgmatt1(1,1),gpizda1(1,1)) + call matmat2(a_temp(1,1),auxgmatt2(1,1),gpizda2(1,1)) + if (shield_mode.eq.0) then fac_shield(i)=1.0d0 fac_shield(j)=1.0d0 @@ -4577,6 +4824,18 @@ if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn3',i,j,0.5d0*(pizda(1,1)+pizda(2,2)) +!C#ifdef NEWCORR +!C Derivatives in theta + gloc(nphi+i,icg)=gloc(nphi+i,icg) & + +0.5d0*(gpizda1(1,1)+gpizda1(2,2))*wturn3& + *fac_shield(i)*fac_shield(j) + gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)& + +0.5d0*(gpizda2(1,1)+gpizda2(2,2))*wturn3& + *fac_shield(i)*fac_shield(j) +!C#endif + + + if ((fac_shield(i).gt.0).and.(fac_shield(j).gt.0).and. & (shield_mode.gt.0)) then !C print *,i,j @@ -4714,8 +4973,14 @@ ! include 'COMMON.CONTROL' real(kind=8),dimension(3) :: ggg real(kind=8),dimension(2,2) :: auxmat,auxmat1,auxmat2,pizda,& - e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2 - real(kind=8),dimension(2) :: auxvec,auxvec1 + e1t,e2t,e3t,e1tder,e2tder,e3tder,e1a,ae3,ae3e2,& + gte1t,gte2t,gte3t,& + gte1a,gtae3,gtae3e2, ae3gte2,& + gtEpizda1,gtEpizda2,gtEpizda3 + + real(kind=8),dimension(2) :: auxvec,auxvec1,auxgEvec1,auxgEvec2,& + auxgEvec3,auxgvec + !el real(kind=8),dimension(3,4) :: agg,aggi,aggi1,aggj,aggj1 real(kind=8),dimension(2,2) :: auxmat3 !el a_temp !el real(kind=8) :: a22,a23,a32,a33,dxi,dyi,dzi,dx_normi,dy_normi,& @@ -4727,7 +4992,7 @@ !el local variables integer :: i,j,iti1,iti2,iti3,l,k,ilist,iresshield real(kind=8) :: eello_turn4,s1,s2,s3,zj,fracinbuf,eello_t4,& - rlocshield + rlocshield,gs23,gs32,gsE13,gs13,gs21,gsE31,gsEE1,gsEE2,gsEE3 j=i+3 ! if (j.ne.20) return @@ -4782,15 +5047,56 @@ call transpose2(EUg(1,1,i+1),e1t(1,1)) call transpose2(Eug(1,1,i+2),e2t(1,1)) call transpose2(Eug(1,1,i+3),e3t(1,1)) +!C Ematrix derivative in theta + call transpose2(gtEUg(1,1,i+1),gte1t(1,1)) + call transpose2(gtEug(1,1,i+2),gte2t(1,1)) + call transpose2(gtEug(1,1,i+3),gte3t(1,1)) + call matmat2(e1t(1,1),a_temp(1,1),e1a(1,1)) call matvec2(e1a(1,1),Ub2(1,i+3),auxvec(1)) + call matmat2(gte1t(1,1),a_temp(1,1),gte1a(1,1)) + call matvec2(e1a(1,1),gUb2(1,i+3),auxgvec(1)) +!c auxalary matrix of E i+1 + call matvec2(gte1a(1,1),Ub2(1,i+3),auxgEvec1(1)) s1=scalar2(b1(1,iti2),auxvec(1)) +!c derivative of theta i+2 with constant i+3 + gs23=scalar2(gtb1(1,i+2),auxvec(1)) +!c derivative of theta i+2 with constant i+2 + gs32=scalar2(b1(1,i+2),auxgvec(1)) +!c derivative of E matix in theta of i+1 + gsE13=scalar2(b1(1,i+2),auxgEvec1(1)) + call matmat2(a_temp(1,1),e3t(1,1),ae3(1,1)) + call matmat2(a_temp(1,1),gte3t(1,1),gtae3(1,1)) call matvec2(ae3(1,1),Ub2(1,i+2),auxvec(1)) +!c auxilary matrix auxgvec of Ub2 with constant E matirx + call matvec2(ae3(1,1),gUb2(1,i+2),auxgvec(1)) +!c auxilary matrix auxgEvec1 of E matix with Ub2 constant + call matvec2(gtae3(1,1),Ub2(1,i+2),auxgEvec3(1)) s2=scalar2(b1(1,iti1),auxvec(1)) +!c derivative of theta i+1 with constant i+3 + gs13=scalar2(gtb1(1,i+1),auxvec(1)) +!c derivative of theta i+2 with constant i+1 + gs21=scalar2(b1(1,i+1),auxgvec(1)) +!c derivative of theta i+3 with constant i+1 + gsE31=scalar2(b1(1,i+1),auxgEvec3(1)) + call matmat2(ae3(1,1),e2t(1,1),ae3e2(1,1)) + call matmat2(gtae3(1,1),e2t(1,1),gtae3e2(1,1)) +!c ae3gte2 is derivative over i+2 + call matmat2(ae3(1,1),gte2t(1,1),ae3gte2(1,1)) + call matmat2(ae3e2(1,1),e1t(1,1),pizda(1,1)) + call matmat2(ae3e2(1,1),gte1t(1,1),gtEpizda1(1,1)) +!c i+2 + call matmat2(ae3gte2(1,1),e1t(1,1),gtEpizda2(1,1)) +!c i+3 + call matmat2(gtae3e2(1,1),e1t(1,1),gtEpizda3(1,1)) + s3=0.5d0*(pizda(1,1)+pizda(2,2)) + gsEE1=0.5d0*(gtEpizda1(1,1)+gtEpizda1(2,2)) + gsEE2=0.5d0*(gtEpizda2(1,1)+gtEpizda2(2,2)) + gsEE3=0.5d0*(gtEpizda3(1,1)+gtEpizda3(2,2)) if (shield_mode.eq.0) then fac_shield(i)=1.0 fac_shield(j)=1.0 @@ -4844,7 +5150,21 @@ ! print *,"gshieldc_t4(k,j+1)",j,gshieldc_t4(k,j+1) enddo endif - +#ifdef NEWCORR + gloc(nphi+i,icg)=gloc(nphi+i,icg)& + -(gs13+gsE13+gsEE1)*wturn4& + *fac_shield(i)*fac_shield(j) + gloc(nphi+i+1,icg)= gloc(nphi+i+1,icg)& + -(gs23+gs21+gsEE2)*wturn4& + *fac_shield(i)*fac_shield(j) + + gloc(nphi+i+2,icg)= gloc(nphi+i+2,icg)& + -(gs32+gsE31+gsEE3)*wturn4& + *fac_shield(i)*fac_shield(j) + +!c gloc(nphi+i+1,icg)=gloc(nphi+i+1,icg)- +!c & gs2 +#endif if (energy_dec) write (iout,'(a6,2i5,0pf7.3)') & 'eturn4',i,j,-(s1+s2+s3) !d write (2,*) 'i,',i,' j',j,'eello_turn4',-(s1+s2+s3), @@ -5897,7 +6217,7 @@ end subroutine theteng #else !----------------------------------------------------------------------------- - subroutine ebend(etheta,ethetacnstr) + subroutine ebend(etheta) ! ! Evaluate the virtual-bond-angle energy given the virtual-bond dihedral ! angles gamma and its derivatives in consecutive thetas and gammas. @@ -6103,34 +6423,6 @@ enddo !-----------thete constrains ! if (tor_mode.ne.2) then - ethetacnstr=0.0d0 -! print *,ithetaconstr_start,ithetaconstr_end,"TU" - do i=ithetaconstr_start,ithetaconstr_end - itheta=itheta_constr(i) - thetiii=theta(itheta) - difi=pinorm(thetiii-theta_constr0(i)) - if (difi.gt.theta_drange(i)) then - difi=difi-theta_drange(i) - ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 - gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & - +for_thet_constr(i)*difi**3 - else if (difi.lt.-drange(i)) then - difi=difi+drange(i) - ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 - gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & - +for_thet_constr(i)*difi**3 - else - difi=0.0 - endif - if (energy_dec) then - write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", & - i,itheta,rad2deg*thetiii, & - rad2deg*theta_constr0(i), rad2deg*theta_drange(i), & - rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, & - gloc(itheta+nphi-2,icg) - endif - enddo -! endif return end subroutine ebend @@ -7019,7 +7311,7 @@ end subroutine etor_d #else !----------------------------------------------------------------------------- - subroutine etor(etors,edihcnstr) + subroutine etor(etors) ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' ! include 'COMMON.VAR' @@ -7101,8 +7393,156 @@ ! write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) enddo ! 6/20/98 - dihedral angle constraints - edihcnstr=0.0d0 -! do i=1,ndih_constr + return + end subroutine etor +!C The rigorous attempt to derive energy function +!------------------------------------------------------------------------------------------- + subroutine etor_kcc(etors) + double precision c1(0:maxval_kcc),c2(0:maxval_kcc) + real(kind=8) :: etors,glocig,glocit1,glocit2,sinthet1,& + sinthet2,costhet1,costhet2,sint1t2,sint1t2n,phii,sinphi,cosphi,& + sint1t2n1,sumvalc,gradvalct1,gradvalct2,sumvals,gradvalst1,& + gradvalst2,etori + logical lprn + integer :: i,j,itori,itori1,nval,k,l + + if (lprn) write (iout,*) "etor_kcc tor_mode",tor_mode + etors=0.0D0 + do i=iphi_start,iphi_end +!C ANY TWO ARE DUMMY ATOMS in row CYCLE +!c if (((itype(i-3).eq.ntyp1).and.(itype(i-2).eq.ntyp1)).or. +!c & ((itype(i-2).eq.ntyp1).and.(itype(i-1).eq.ntyp1)) .or. +!c & ((itype(i-1).eq.ntyp1).and.(itype(i).eq.ntyp1))) cycle + if (itype(i-2,1).eq.ntyp1.or. itype(i-1,1).eq.ntyp1 & + .or. itype(i,1).eq.ntyp1 .or. itype(i-3,1).eq.ntyp1) cycle + itori=itortyp(itype(i-2,1)) + itori1=itortyp(itype(i-1,1)) + phii=phi(i) + glocig=0.0D0 + glocit1=0.0d0 + glocit2=0.0d0 +!C to avoid multiple devision by 2 +!c theti22=0.5d0*theta(i) +!C theta 12 is the theta_1 /2 +!C theta 22 is theta_2 /2 +!c theti12=0.5d0*theta(i-1) +!C and appropriate sinus function + sinthet1=dsin(theta(i-1)) + sinthet2=dsin(theta(i)) + costhet1=dcos(theta(i-1)) + costhet2=dcos(theta(i)) +!C to speed up lets store its mutliplication + sint1t2=sinthet2*sinthet1 + sint1t2n=1.0d0 +!C \sum_{i=1}^n (sin(theta_1) * sin(theta_2))^n * (c_n* cos(n*gamma) +!C +d_n*sin(n*gamma)) * +!C \sum_{i=1}^m (1+a_m*Tb_m(cos(theta_1 /2))+b_m*Tb_m(cos(theta_2 /2))) +!C we have two sum 1) Non-Chebyshev which is with n and gamma + nval=nterm_kcc_Tb(itori,itori1) + c1(0)=0.0d0 + c2(0)=0.0d0 + c1(1)=1.0d0 + c2(1)=1.0d0 + do j=2,nval + c1(j)=c1(j-1)*costhet1 + c2(j)=c2(j-1)*costhet2 + enddo + etori=0.0d0 + + do j=1,nterm_kcc(itori,itori1) + cosphi=dcos(j*phii) + sinphi=dsin(j*phii) + sint1t2n1=sint1t2n + sint1t2n=sint1t2n*sint1t2 + sumvalc=0.0d0 + gradvalct1=0.0d0 + gradvalct2=0.0d0 + do k=1,nval + do l=1,nval + sumvalc=sumvalc+v1_kcc(l,k,j,itori1,itori)*c1(k)*c2(l) + gradvalct1=gradvalct1+ & + (k-1)*v1_kcc(l,k,j,itori1,itori)*c1(k-1)*c2(l) + gradvalct2=gradvalct2+ & + (l-1)*v1_kcc(l,k,j,itori1,itori)*c1(k)*c2(l-1) + enddo + enddo + gradvalct1=-gradvalct1*sinthet1 + gradvalct2=-gradvalct2*sinthet2 + sumvals=0.0d0 + gradvalst1=0.0d0 + gradvalst2=0.0d0 + do k=1,nval + do l=1,nval + sumvals=sumvals+v2_kcc(l,k,j,itori1,itori)*c1(k)*c2(l) + gradvalst1=gradvalst1+ & + (k-1)*v2_kcc(l,k,j,itori1,itori)*c1(k-1)*c2(l) + gradvalst2=gradvalst2+ & + (l-1)*v2_kcc(l,k,j,itori1,itori)*c1(k)*c2(l-1) + enddo + enddo + gradvalst1=-gradvalst1*sinthet1 + gradvalst2=-gradvalst2*sinthet2 + if (lprn) write (iout,*)j,"sumvalc",sumvalc," sumvals",sumvals + etori=etori+sint1t2n*(sumvalc*cosphi+sumvals*sinphi) +!C glocig is the gradient local i site in gamma + glocig=glocig+j*sint1t2n*(sumvals*cosphi-sumvalc*sinphi) +!C now gradient over theta_1 + glocit1=glocit1+sint1t2n*(gradvalct1*cosphi+gradvalst1*sinphi)& + +j*sint1t2n1*costhet1*sinthet2*(sumvalc*cosphi+sumvals*sinphi) + glocit2=glocit2+sint1t2n*(gradvalct2*cosphi+gradvalst2*sinphi)& + +j*sint1t2n1*sinthet1*costhet2*(sumvalc*cosphi+sumvals*sinphi) + enddo ! j + etors=etors+etori + gloc(i-3,icg)=gloc(i-3,icg)+wtor*glocig +!C derivative over theta1 + gloc(nphi+i-3,icg)=gloc(nphi+i-3,icg)+wtor*glocit1 +!C now derivative over theta2 + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wtor*glocit2 + if (lprn) then + write (iout,*) i-2,i-1,itype(i-2,1),itype(i-1,1),itori,itori1,& + theta(i-1)*rad2deg,theta(i)*rad2deg,phii*rad2deg,etori + write (iout,*) "c1",(c1(k),k=0,nval), & + " c2",(c2(k),k=0,nval) + endif + enddo + return + end subroutine etor_kcc +!------------------------------------------------------------------------------ + + subroutine etor_constr(edihcnstr) + real(kind=8) :: etors,edihcnstr + logical :: lprn +!el local variables + integer :: i,j,iblock,itori,itori1 + real(kind=8) :: phii,gloci,v1ij,v2ij,cosphi,sinphi,& + vl1ij,vl2ij,vl3ij,pom1,difi,etors_ii,pom,& + gaudih_i,gauder_i,s,cos_i,dexpcos_i + + if (raw_psipred) then + do i=idihconstr_start,idihconstr_end + itori=idih_constr(i) + phii=phi(itori) + gaudih_i=vpsipred(1,i) + gauder_i=0.0d0 + do j=1,2 + s = sdihed(j,i) + cos_i=(1.0d0-dcos(phii-phibound(j,i)))/s**2 + dexpcos_i=dexp(-cos_i*cos_i) + gaudih_i=gaudih_i+vpsipred(j+1,i)*dexpcos_i + gauder_i=gauder_i-2*vpsipred(j+1,i)*dsin(phii-phibound(j,i)) & + *cos_i*dexpcos_i/s**2 + enddo + edihcnstr=edihcnstr-wdihc*dlog(gaudih_i) + gloc(itori-3,icg)=gloc(itori-3,icg)-wdihc*gauder_i/gaudih_i + if (energy_dec) & + write (iout,'(2i5,f8.3,f8.5,2(f8.5,2f8.3),f10.5)') & + i,itori,phii*rad2deg,vpsipred(1,i),vpsipred(2,i),& + phibound(1,i)*rad2deg,sdihed(1,i)*rad2deg,vpsipred(3,i),& + phibound(2,i)*rad2deg,sdihed(2,i)*rad2deg,& + -wdihc*dlog(gaudih_i) + enddo + else + do i=idihconstr_start,idihconstr_end itori=idih_constr(i) phii=phi(itori) @@ -7118,13 +7558,13 @@ else difi=0.0 endif -!d write (iout,'(2i5,4f8.3,2e14.5)') i,itori,rad2deg*phii, -!d & rad2deg*phi0(i), rad2deg*drange(i), -!d & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) enddo -!d write (iout,*) 'edihcnstr',edihcnstr + + endif + return - end subroutine etor + + end subroutine etor_constr !----------------------------------------------------------------------------- subroutine etor_d(etors_d) ! 6/23/01 Compute double torsional energy @@ -7215,6 +7655,75 @@ return end subroutine etor_d #endif + + subroutine ebend_kcc(etheta) + logical lprn + double precision thybt1(maxang_kcc),etheta + integer :: i,iti,j,ihelp + real (kind=8) :: sinthet,costhet,sumth1thyb,gradthybt1 +!C Set lprn=.true. for debugging + lprn=energy_dec +!c lprn=.true. +!C print *,"wchodze kcc" + if (lprn) write (iout,*) "ebend_kcc tor_mode",tor_mode + etheta=0.0D0 + do i=ithet_start,ithet_end +!c print *,i,itype(i-1),itype(i),itype(i-2) + if ((itype(i-1,1).eq.ntyp1).or.itype(i-2,1).eq.ntyp1 & + .or.itype(i,1).eq.ntyp1) cycle + iti=iabs(itortyp(itype(i-1,1))) + sinthet=dsin(theta(i)) + costhet=dcos(theta(i)) + do j=1,nbend_kcc_Tb(iti) + thybt1(j)=v1bend_chyb(j,iti) + enddo + sumth1thyb=v1bend_chyb(0,iti)+ & + tschebyshev(1,nbend_kcc_Tb(iti),thybt1(1),costhet) + if (lprn) write (iout,*) i-1,itype(i-1,1),iti,theta(i)*rad2deg,& + sumth1thyb + ihelp=nbend_kcc_Tb(iti)-1 + gradthybt1=gradtschebyshev(0,ihelp,thybt1(1),costhet) + etheta=etheta+sumth1thyb +!C print *,sumth1thyb,gradthybt1,sinthet*(-0.5d0) + gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)-wang*gradthybt1*sinthet + enddo + return + end subroutine ebend_kcc +!c------------ +!c------------------------------------------------------------------------------------- + subroutine etheta_constr(ethetacnstr) + real (kind=8) :: ethetacnstr,thetiii,difi + integer :: i,itheta + ethetacnstr=0.0d0 +!C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=ithetaconstr_start,ithetaconstr_end + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & + +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) & + +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc",& + i,itheta,rad2deg*thetiii,& + rad2deg*theta_constr0(i), rad2deg*theta_drange(i),& + rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4,& + gloc(itheta+nphi-2,icg) + endif + enddo + return + end subroutine etheta_constr + !----------------------------------------------------------------------------- subroutine eback_sc_corr(esccor) ! 7/21/2007 Correlations between the backbone-local and side-chain-local @@ -8411,9 +8920,9 @@ iti1 = itortyp(itype(i+1,1)) if (j.lt.nres-1) then - itj1 = itortyp(itype(j+1,1)) + itj1 = itype2loc(itype(j+1,1)) else - itj1=ntortyp+1 + itj1=nloctyp endif do iii=1,2 dipi(iii,1)=Ub2(iii,i) @@ -16123,20 +16632,66 @@ ! ! Calculate the virtual-bond-angle energy. ! - call ebend(ebe,ethetacnstr) -! ! Calculate the SC local energy. ! call vec_and_deriv call esc(escloc) ! + if (wang.gt.0d0) then + if (tor_mode.eq.0) then + call ebend(ebe) + else +!C ebend kcc is Kubo cumulant clustered rigorous attemp to derive the +!C energy function + call ebend_kcc(ebe) + endif + else + ebe=0.0d0 + endif + ethetacnstr=0.0d0 + if (with_theta_constr) call etheta_constr(ethetacnstr) + +! write(iout,*) "in etotal afer ebe",ipot + +! print *,"Processor",myrank," computed UB" +! +! Calculate the SC local energy. +! + call esc(escloc) +!elwrite(iout,*) "in etotal afer esc",ipot +! print *,"Processor",myrank," computed USC" +! +! Calculate the virtual-bond torsional energy. +! +!d print *,'nterm=',nterm +! if (wtor.gt.0) then +! call etor(etors,edihcnstr) +! else +! etors=0 +! edihcnstr=0 +! endif + if (wtor.gt.0.0d0) then + if (tor_mode.eq.0) then + call etor(etors) + else +!C etor kcc is Kubo cumulant clustered rigorous attemp to derive the +!C energy function + call etor_kcc(etors) + endif + else + etors=0.0d0 + endif + edihcnstr=0.0d0 + if (ndih_constr.gt.0) call etor_constr(edihcnstr) + ! Calculate the virtual-bond torsional energy. ! - call etor(etors,edihcnstr) ! ! 6/23/01 Calculate double-torsional energy ! + if ((wtor_d.gt.0.0d0).and.(tor_mode.eq.0)) then call etor_d(etors_d) + endif ! ! 21/5/07 Calculate local sicdechain correlation energy ! @@ -19958,6 +20513,24 @@ allocate(Ug2DtEUgder(2,2,2,nres)) allocate(DtUg2EUgder(2,2,2,nres)) !(2,2,2,maxres) + allocate(b1(2,nres)) !(2,-maxtor:maxtor) + allocate(b2(2,nres)) !(2,-maxtor:maxtor) + allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor) + allocate(b2tilde(2,nres)) !(2,-maxtor:maxtor) + + allocate(ctilde(2,2,nres)) + allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor) + allocate(gtb1(2,nres)) + allocate(gtb2(2,nres)) + allocate(cc(2,2,nres)) + allocate(dd(2,2,nres)) + allocate(ee(2,2,nres)) + allocate(gtcc(2,2,nres)) + allocate(gtdd(2,2,nres)) + allocate(gtee(2,2,nres)) + allocate(gUb2(2,nres)) + allocate(gteUg(2,2,nres)) + ! common /rotat_old/ allocate(costab(nres)) allocate(sintab(nres)) @@ -25974,4 +26547,50 @@ dPOLdOM2 = 0.0d0 RETURN END SUBROUTINE elgrad_init + + double precision function tschebyshev(m,n,x,y) + implicit none + integer i,m,n + double precision x(n),y,yy(0:maxvar),aux +!c Tschebyshev polynomial. Note that the first term is omitted +!c m=0: the constant term is included +!c m=1: the constant term is not included + yy(0)=1.0d0 + yy(1)=y + do i=2,n + yy(i)=2*yy(1)*yy(i-1)-yy(i-2) + enddo + aux=0.0d0 + do i=m,n + aux=aux+x(i)*yy(i) + enddo + tschebyshev=aux + return + end function tschebyshev +!C-------------------------------------------------------------------------- + double precision function gradtschebyshev(m,n,x,y) + implicit none + integer i,m,n + double precision x(n+1),y,yy(0:maxvar),aux +!c Tschebyshev polynomial. Note that the first term is omitted +!c m=0: the constant term is included +!c m=1: the constant term is not included + yy(0)=1.0d0 + yy(1)=2.0d0*y + do i=2,n + yy(i)=2*y*yy(i-1)-yy(i-2) + enddo + aux=0.0d0 + do i=m,n + aux=aux+x(i+1)*yy(i)*(i+1) +!C print *, x(i+1),yy(i),i + enddo + gradtschebyshev=aux + return + end function gradtschebyshev + + + + + end module energy diff --git a/source/unres/io.F90 b/source/unres/io.F90 index cc8b681..abce524 100644 --- a/source/unres/io.F90 +++ b/source/unres/io.F90 @@ -736,10 +736,12 @@ real(kind=8),dimension(3,maxres2+2) :: c_alloc real(kind=8),dimension(3,0:maxres2) :: dc_alloc + real(kind=8),dimension(:,:), allocatable :: secprob integer,dimension(maxres) :: itype_alloc integer :: iti,nsi,maxsi,itrial,itmp - real(kind=8) :: wlong,scalscp,co,ssscale + real(kind=8) :: wlong,scalscp,co,ssscale,phihel,phibet,sigmahel,& + sigmabet,sumv allocate(weights(n_ene)) !----------------------------- allocate(c(3,2*maxres+2)) !(3,maxres2+2) maxres2=2*maxres @@ -1155,6 +1157,7 @@ enddo read (inp,*) ndih_constr if (ndih_constr.gt.0) then + raw_psipred=.false. allocate(idih_constr(ndih_constr),idih_nconstr(ndih_constr)) !(maxdih_constr) allocate(phi0(ndih_constr),drange(ndih_constr)) !(maxdih_constr) allocate(ftors(ndih_constr)) !(maxdih_constr) @@ -1181,6 +1184,45 @@ phibound(1,ii) = phi0(i)-drange(i) phibound(2,ii) = phi0(i)+drange(i) enddo + else if (ndih_constr.lt.0) then + raw_psipred=.true. + allocate(secprob(3,nres)) + allocate(vpsipred(3,nres)) + allocate(sdihed(2,nres)) + call card_concat(weightcard,.true.) + call reada(weightcard,"PHIHEL",phihel,50.0D0) + call reada(weightcard,"PHIBET",phibet,180.0D0) + call reada(weightcard,"SIGMAHEL",sigmahel,30.0d0) + call reada(weightcard,"SIGMABET",sigmabet,40.0d0) + call reada(weightcard,"WDIHC",wdihc,0.591D0) + write (iout,*) "Weight of dihedral angle restraints",wdihc + read(inp,'(9x,3f7.3)') & + (secprob(1,i),secprob(2,i),secprob(3,i),i=nnt,nct) + write (iout,*) "The secprob array" + do i=nnt,nct + write (iout,'(i5,3f8.3)') i,(secprob(j,i),j=1,3) + enddo + ndih_constr=0 + do i=nnt+3,nct + if (itype(i-3,1).ne.ntyp1 .and. itype(i-2,1).ne.ntyp1 & + .and. itype(i-1,1).ne.ntyp1 .and. itype(i,1).ne.ntyp1) then + ndih_constr=ndih_constr+1 + idih_constr(ndih_constr)=i + sumv=0.0d0 + do j=1,3 + vpsipred(j,ndih_constr)=secprob(j,i-1)*secprob(j,i-2) + sumv=sumv+vpsipred(j,ndih_constr) + enddo + do j=1,3 + vpsipred(j,ndih_constr)=vpsipred(j,ndih_constr)/sumv + enddo + phibound(1,ndih_constr)=phihel*deg2rad + phibound(2,ndih_constr)=phibet*deg2rad + sdihed(1,ndih_constr)=sigmahel*deg2rad + sdihed(2,ndih_constr)=sigmabet*deg2rad + endif + enddo + endif if (with_theta_constr) then !C with_theta_constr is keyword allowing for occurance of theta constrains @@ -1371,7 +1413,10 @@ enddo return else - call read_angles(inp,*36) + write(iout,*) "read angles from input" + call read_angles(inp,*36) + call chainbuild + endif goto 37 36 write (iout,'(a)') 'Error reading angle file.' diff --git a/source/unres/io_config.F90 b/source/unres/io_config.F90 index 2b1db84..a31dead 100644 --- a/source/unres/io_config.F90 +++ b/source/unres/io_config.F90 @@ -734,19 +734,21 @@ character(len=1) :: t1,t2,t3 character(len=1) :: onelett(4) = (/"G","A","P","D"/) character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/) - logical :: lprint,LaTeX + logical :: lprint,LaTeX,SPLIT_FOURIERTOR real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob) real(kind=8),dimension(13) :: buse character(len=3) :: lancuch !,ucase !el local variables - integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm - integer :: maxinter,junk,kk,ii,ncatprotparm + integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj + integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,& dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,& sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,& res1,epsijlip,epspeptube,epssctube,sigmapeptube, & sigmasctube - integer :: ichir1,ichir2 + integer :: ichir1,ichir2,ijunk + character*3 string + ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)) @@ -1017,6 +1019,7 @@ ! Read the parameters of Utheta determined from ab initio surfaces ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203 ! + IF (tor_mode.eq.0) THEN read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,& ntheterm3,nsingle,ndouble nntheterm=max0(ntheterm,ntheterm2,ntheterm3) @@ -1209,6 +1212,29 @@ enddo call flush(iout) endif + ELSE +!C here will be the apropriate recalibrating for D-aminoacid + read (ithep,*,end=121,err=121) nthetyp + allocate(nbend_kcc_Tb(-nthetyp:nthetyp)) + allocate(v1bend_chyb(0:36,-nthetyp:nthetyp)) + do i=-nthetyp+1,nthetyp-1 + read (ithep,*,end=121,err=121) nbend_kcc_Tb(i) + do j=0,nbend_kcc_Tb(i) + read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i) + enddo + enddo + if (lprint) then + write (iout,'(a)') & + "Parameters of the valence-only potentials" + do i=-nthetyp+1,nthetyp-1 + write (iout,'(2a)') "Type ",toronelet(i) + do k=0,nbend_kcc_Tb(i) + write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i) + enddo + enddo + endif + ENDIF ! TOR_MODE + write (2,*) "Start reading THETA_PDB",ithep_pdb do i=1,ntyp ! write (2,*) 'i=',i @@ -1501,6 +1527,442 @@ #endif close(irotam) + +!C +!C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local +!C interaction energy of the Gly, Ala, and Pro prototypes. +!C + read (ifourier,*) nloctyp + SPLIT_FOURIERTOR = nloctyp.lt.0 + nloctyp = iabs(nloctyp) +!C allocate(b1(2,nres)) !(2,-maxtor:maxtor) +!C allocate(b2(2,nres)) !(2,-maxtor:maxtor) +!C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor) +!C allocate(ctilde(2,2,nres)) +!C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor) +!C allocate(gtb1(2,nres)) +!C allocate(gtb2(2,nres)) +!C allocate(cc(2,2,nres)) +!C allocate(dd(2,2,nres)) +!C allocate(ee(2,2,nres)) +!C allocate(gtcc(2,2,nres)) +!C allocate(gtdd(2,2,nres)) +!C allocate(gtee(2,2,nres)) + +#ifdef NEWCORR + allocate(itype2loc(-ntyp1:ntyp1)) + allocate(iloctyp(-nloctyp:nloctyp)) + allocate(bnew1(3,2,-nloctyp:nloctyp)) + allocate(bnew2(3,2,-nloctyp:nloctyp)) + allocate(ccnew(3,2,-nloctyp:nloctyp)) + allocate(ddnew(3,2,-nloctyp:nloctyp)) + allocate(e0new(3,-nloctyp:nloctyp)) + allocate(eenew(2,2,2,-nloctyp:nloctyp)) + allocate(bnew1tor(3,2,-nloctyp:nloctyp)) + allocate(bnew2tor(3,2,-nloctyp:nloctyp)) + allocate(ccnewtor(3,2,-nloctyp:nloctyp)) + allocate(ddnewtor(3,2,-nloctyp:nloctyp)) + allocate(e0newtor(3,-nloctyp:nloctyp)) + allocate(eenewtor(2,2,2,-nloctyp:nloctyp)) + + read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp) + read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1) + itype2loc(ntyp1)=nloctyp + iloctyp(nloctyp)=ntyp1 + do i=1,ntyp1 + itype2loc(-i)=-itype2loc(i) + enddo +#else + iloctyp(0)=10 + iloctyp(1)=9 + iloctyp(2)=20 + iloctyp(3)=ntyp1 +#endif + do i=1,nloctyp + iloctyp(-i)=-iloctyp(i) + enddo +!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) +!c write (iout,*) "nloctyp",nloctyp, +!c & " iloctyp",(iloctyp(i),i=0,nloctyp) +!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1) +!c write (iout,*) "nloctyp",nloctyp, +!c & " iloctyp",(iloctyp(i),i=0,nloctyp) +#ifdef NEWCORR + do i=0,nloctyp-1 +!c write (iout,*) "NEWCORR",i + read (ifourier,*,end=115,err=115) + do ii=1,3 + do j=1,2 + read (ifourier,*,end=115,err=115) bnew1(ii,j,i) + enddo + enddo +!c write (iout,*) "NEWCORR BNEW1" +!c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2) + do ii=1,3 + do j=1,2 + read (ifourier,*,end=115,err=115) bnew2(ii,j,i) + enddo + enddo +!c write (iout,*) "NEWCORR BNEW2" +!c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2) + do kk=1,3 + read (ifourier,*,end=115,err=115) ccnew(kk,1,i) + read (ifourier,*,end=115,err=115) ccnew(kk,2,i) + enddo +!c write (iout,*) "NEWCORR CCNEW" +!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2) + do kk=1,3 + read (ifourier,*,end=115,err=115) ddnew(kk,1,i) + read (ifourier,*,end=115,err=115) ddnew(kk,2,i) + enddo +!c write (iout,*) "NEWCORR DDNEW" +!c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2) + do ii=1,2 + do jj=1,2 + do kk=1,2 + read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i) + enddo + enddo + enddo +!c write (iout,*) "NEWCORR EENEW1" +!c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2) + do ii=1,3 + read (ifourier,*,end=115,err=115) e0new(ii,i) + enddo +!c write (iout,*) (e0new(ii,i),ii=1,3) + enddo +!c write (iout,*) "NEWCORR EENEW" + do i=0,nloctyp-1 + do ii=1,3 + ccnew(ii,1,i)=ccnew(ii,1,i)/2 + ccnew(ii,2,i)=ccnew(ii,2,i)/2 + ddnew(ii,1,i)=ddnew(ii,1,i)/2 + ddnew(ii,2,i)=ddnew(ii,2,i)/2 + enddo + enddo + do i=1,nloctyp-1 + do ii=1,3 + bnew1(ii,1,-i)= bnew1(ii,1,i) + bnew1(ii,2,-i)=-bnew1(ii,2,i) + bnew2(ii,1,-i)= bnew2(ii,1,i) + bnew2(ii,2,-i)=-bnew2(ii,2,i) + enddo + do ii=1,3 +!c ccnew(ii,1,i)=ccnew(ii,1,i)/2 +!c ccnew(ii,2,i)=ccnew(ii,2,i)/2 +!c ddnew(ii,1,i)=ddnew(ii,1,i)/2 +!c ddnew(ii,2,i)=ddnew(ii,2,i)/2 + ccnew(ii,1,-i)=ccnew(ii,1,i) + ccnew(ii,2,-i)=-ccnew(ii,2,i) + ddnew(ii,1,-i)=ddnew(ii,1,i) + ddnew(ii,2,-i)=-ddnew(ii,2,i) + enddo + e0new(1,-i)= e0new(1,i) + e0new(2,-i)=-e0new(2,i) + e0new(3,-i)=-e0new(3,i) + do kk=1,2 + eenew(kk,1,1,-i)= eenew(kk,1,1,i) + eenew(kk,1,2,-i)=-eenew(kk,1,2,i) + eenew(kk,2,1,-i)=-eenew(kk,2,1,i) + eenew(kk,2,2,-i)= eenew(kk,2,2,i) + enddo + enddo + if (lprint) then + write (iout,'(a)') "Coefficients of the multibody terms" + do i=-nloctyp+1,nloctyp-1 + write (iout,*) "Type: ",onelet(iloctyp(i)) + write (iout,*) "Coefficients of the expansion of B1" + do j=1,2 + write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3) + enddo + write (iout,*) "Coefficients of the expansion of B2" + do j=1,2 + write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3) + enddo + write (iout,*) "Coefficients of the expansion of C" + write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3) + write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3) + write (iout,*) "Coefficients of the expansion of D" + write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3) + write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3) + write (iout,*) "Coefficients of the expansion of E" + write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3) + do j=1,2 + do k=1,2 + write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2) + enddo + enddo + enddo + endif + IF (SPLIT_FOURIERTOR) THEN + do i=0,nloctyp-1 +!c write (iout,*) "NEWCORR TOR",i + read (ifourier,*,end=115,err=115) + do ii=1,3 + do j=1,2 + read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i) + enddo + enddo +!c write (iout,*) "NEWCORR BNEW1 TOR" +!c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2) + do ii=1,3 + do j=1,2 + read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i) + enddo + enddo +!c write (iout,*) "NEWCORR BNEW2 TOR" +!c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2) + do kk=1,3 + read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i) + read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i) + enddo +!c write (iout,*) "NEWCORR CCNEW TOR" +!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2) + do kk=1,3 + read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i) + read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i) + enddo +!c write (iout,*) "NEWCORR DDNEW TOR" +!c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2) + do ii=1,2 + do jj=1,2 + do kk=1,2 + read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i) + enddo + enddo + enddo +!c write (iout,*) "NEWCORR EENEW1 TOR" +!c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2) + do ii=1,3 + read (ifourier,*,end=115,err=115) e0newtor(ii,i) + enddo +!c write (iout,*) (e0newtor(ii,i),ii=1,3) + enddo +!c write (iout,*) "NEWCORR EENEW TOR" + do i=0,nloctyp-1 + do ii=1,3 + ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2 + ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2 + ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2 + ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2 + enddo + enddo + do i=1,nloctyp-1 + do ii=1,3 + bnew1tor(ii,1,-i)= bnew1tor(ii,1,i) + bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i) + bnew2tor(ii,1,-i)= bnew2tor(ii,1,i) + bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i) + enddo + do ii=1,3 + ccnewtor(ii,1,-i)=ccnewtor(ii,1,i) + ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i) + ddnewtor(ii,1,-i)=ddnewtor(ii,1,i) + ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i) + enddo + e0newtor(1,-i)= e0newtor(1,i) + e0newtor(2,-i)=-e0newtor(2,i) + e0newtor(3,-i)=-e0newtor(3,i) + do kk=1,2 + eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i) + eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i) + eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i) + eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i) + enddo + enddo + if (lprint) then + write (iout,'(a)') & + "Single-body coefficients of the torsional potentials" + do i=-nloctyp+1,nloctyp-1 + write (iout,*) "Type: ",onelet(iloctyp(i)) + write (iout,*) "Coefficients of the expansion of B1tor" + do j=1,2 + write (iout,'(3hB1(,i1,1h),3f10.5)') & + j,(bnew1tor(k,j,i),k=1,3) + enddo + write (iout,*) "Coefficients of the expansion of B2tor" + do j=1,2 + write (iout,'(3hB2(,i1,1h),3f10.5)') & + j,(bnew2tor(k,j,i),k=1,3) + enddo + write (iout,*) "Coefficients of the expansion of Ctor" + write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3) + write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3) + write (iout,*) "Coefficients of the expansion of Dtor" + write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3) + write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3) + write (iout,*) "Coefficients of the expansion of Etor" + write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3) + do j=1,2 + do k=1,2 + write (iout,'(1hE,2i1,2f10.5)') & + j,k,(eenewtor(l,j,k,i),l=1,2) + enddo + enddo + enddo + endif + ELSE + do i=-nloctyp+1,nloctyp-1 + do ii=1,3 + do j=1,2 + bnew1tor(ii,j,i)=bnew1(ii,j,i) + enddo + enddo + do ii=1,3 + do j=1,2 + bnew2tor(ii,j,i)=bnew2(ii,j,i) + enddo + enddo + do ii=1,3 + ccnewtor(ii,1,i)=ccnew(ii,1,i) + ccnewtor(ii,2,i)=ccnew(ii,2,i) + ddnewtor(ii,1,i)=ddnew(ii,1,i) + ddnewtor(ii,2,i)=ddnew(ii,2,i) + enddo + enddo + ENDIF !SPLIT_FOURIER_TOR +#else + allocate(ccold(2,2,-nloctyp-1:nloctyp+1)) + allocate(ddold(2,2,-nloctyp-1:nloctyp+1)) + allocate(eeold(2,2,-nloctyp-1:nloctyp+1)) + + if (lprint) & + write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)" + do i=0,nloctyp-1 + read (ifourier,*,end=115,err=115) + read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13) + if (lprint) then + write (iout,*) 'Type ',onelet(iloctyp(i)) + write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13) + endif + if (i.gt.0) then + b(2,-i)= b(2,i) + b(3,-i)= b(3,i) + b(4,-i)=-b(4,i) + b(5,-i)=-b(5,i) + endif +c B1(1,i) = b(3) +c B1(2,i) = b(5) +c B1(1,-i) = b(3) +c B1(2,-i) = -b(5) +c b1(1,i)=0.0d0 +c b1(2,i)=0.0d0 +c B1tilde(1,i) = b(3) +c B1tilde(2,i) =-b(5) +c B1tilde(1,-i) =-b(3) +c B1tilde(2,-i) =b(5) +c b1tilde(1,i)=0.0d0 +c b1tilde(2,i)=0.0d0 +c B2(1,i) = b(2) +c B2(2,i) = b(4) +c B2(1,-i) =b(2) +c B2(2,-i) =-b(4) +cc B1tilde(1,i) = b(3,i) +cc B1tilde(2,i) =-b(5,i) +C B1tilde(1,-i) =-b(3,i) +C B1tilde(2,-i) =b(5,i) +cc b1tilde(1,i)=0.0d0 +cc b1tilde(2,i)=0.0d0 +cc B2(1,i) = b(2,i) +cc B2(2,i) = b(4,i) +C B2(1,-i) =b(2,i) +C B2(2,-i) =-b(4,i) + +c b2(1,i)=0.0d0 +c b2(2,i)=0.0d0 + CCold(1,1,i)= b(7,i) + CCold(2,2,i)=-b(7,i) + CCold(2,1,i)= b(9,i) + CCold(1,2,i)= b(9,i) + CCold(1,1,-i)= b(7,i) + CCold(2,2,-i)=-b(7,i) + CCold(2,1,-i)=-b(9,i) + CCold(1,2,-i)=-b(9,i) +c CC(1,1,i)=0.0d0 +c CC(2,2,i)=0.0d0 +c CC(2,1,i)=0.0d0 +c CC(1,2,i)=0.0d0 +c Ctilde(1,1,i)= CCold(1,1,i) +c Ctilde(1,2,i)= CCold(1,2,i) +c Ctilde(2,1,i)=-CCold(2,1,i) +c Ctilde(2,2,i)=-CCold(2,2,i) +c CC(1,1,i)=0.0d0 +c CC(2,2,i)=0.0d0 +c CC(2,1,i)=0.0d0 +c CC(1,2,i)=0.0d0 +c Ctilde(1,1,i)= CCold(1,1,i) +c Ctilde(1,2,i)= CCold(1,2,i) +c Ctilde(2,1,i)=-CCold(2,1,i) +c Ctilde(2,2,i)=-CCold(2,2,i) + +c Ctilde(1,1,i)=0.0d0 +c Ctilde(1,2,i)=0.0d0 +c Ctilde(2,1,i)=0.0d0 +c Ctilde(2,2,i)=0.0d0 + DDold(1,1,i)= b(6,i) + DDold(2,2,i)=-b(6,i) + DDold(2,1,i)= b(8,i) + DDold(1,2,i)= b(8,i) + DDold(1,1,-i)= b(6,i) + DDold(2,2,-i)=-b(6,i) + DDold(2,1,-i)=-b(8,i) + DDold(1,2,-i)=-b(8,i) +c DD(1,1,i)=0.0d0 +c DD(2,2,i)=0.0d0 +c DD(2,1,i)=0.0d0 +c DD(1,2,i)=0.0d0 +c Dtilde(1,1,i)= DD(1,1,i) +c Dtilde(1,2,i)= DD(1,2,i) +c Dtilde(2,1,i)=-DD(2,1,i) +c Dtilde(2,2,i)=-DD(2,2,i) + +c Dtilde(1,1,i)=0.0d0 +c Dtilde(1,2,i)=0.0d0 +c Dtilde(2,1,i)=0.0d0 +c Dtilde(2,2,i)=0.0d0 + EEold(1,1,i)= b(10,i)+b(11,i) + EEold(2,2,i)=-b(10,i)+b(11,i) + EEold(2,1,i)= b(12,i)-b(13,i) + EEold(1,2,i)= b(12,i)+b(13,i) + EEold(1,1,-i)= b(10,i)+b(11,i) + EEold(2,2,-i)=-b(10,i)+b(11,i) + EEold(2,1,-i)=-b(12,i)+b(13,i) + EEold(1,2,-i)=-b(12,i)-b(13,i) + write(iout,*) "TU DOCHODZE" + print *,"JESTEM" +c ee(1,1,i)=1.0d0 +c ee(2,2,i)=1.0d0 +c ee(2,1,i)=0.0d0 +c ee(1,2,i)=0.0d0 +c ee(2,1,i)=ee(1,2,i) + enddo + if (lprint) then + write (iout,*) + write (iout,*) & + "Coefficients of the cumulants (independent of valence angles)" + do i=-nloctyp+1,nloctyp-1 + write (iout,*) 'Type ',onelet(iloctyp(i)) + write (iout,*) 'B1' + write(iout,'(2f10.5)') B(3,i),B(5,i) + write (iout,*) 'B2' + write(iout,'(2f10.5)') B(2,i),B(4,i) + write (iout,*) 'CC' + do j=1,2 + write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i) + enddo + write(iout,*) 'DD' + do j=1,2 + write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i) + enddo + write(iout,*) 'EE' + do j=1,2 + write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i) + enddo + enddo + endif +#endif + + #ifdef CRYST_TOR ! ! Read torsional parameters in old format @@ -1537,6 +1999,7 @@ ! ! Read torsional parameters ! + IF (TOR_MODE.eq.0) THEN allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) read (itorp,*,end=113,err=113) ntortyp !el from energy module--------- @@ -1722,6 +2185,113 @@ enddo enddo endif + ELSE IF (TOR_MODE.eq.1) THEN + +!C read valence-torsional parameters + read (itorp,*,end=121,err=121) ntortyp + nkcctyp=ntortyp + write (iout,*) "Valence-torsional parameters read in ntortyp",& + ntortyp + read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp) + write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp) +#ifndef NEWCORR + do i=1,ntyp1 + itype2loc(i)=itortyp(i) + enddo +#endif + do i=-ntyp,-1 + itortyp(i)=-itortyp(-i) + enddo + do i=-ntortyp+1,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 +!C first we read the cos and sin gamma parameters + read (itorp,'(13x,a)',end=121,err=121) string + write (iout,*) i,j,string + read (itorp,*,end=121,err=121) & + nterm_kcc(j,i),nterm_kcc_Tb(j,i) +!C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i) + do k=1,nterm_kcc(j,i) + do l=1,nterm_kcc_Tb(j,i) + do ll=1,nterm_kcc_Tb(j,i) + read (itorp,*,end=121,err=121) ii,jj,kk, & + v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) + enddo + enddo + enddo + enddo + enddo + ELSE +#ifdef NEWCORR +!c AL 4/8/16: Calculate coefficients from one-body parameters + ntortyp=nloctyp + allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1) + allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1)) + allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1)) + allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1)) + allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1)) + + do i=-ntyp1,ntyp1 + print *,i,itortyp(i) + itortyp(i)=itype2loc(i) + enddo + write (iout,*) & + "Val-tor parameters calculated from cumulant coefficients ntortyp"& + ,ntortyp + do i=-ntortyp+1,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + nterm_kcc(j,i)=2 + nterm_kcc_Tb(j,i)=3 + do k=1,nterm_kcc_Tb(j,i) + do l=1,nterm_kcc_Tb(j,i) + v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)& + +bnew1tor(k,2,i)*bnew2tor(l,2,j) + v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)& + +bnew1tor(k,2,i)*bnew2tor(l,1,j) + enddo + enddo + do k=1,nterm_kcc_Tb(j,i) + do l=1,nterm_kcc_Tb(j,i) +#ifdef CORRCD + v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) & + -ccnewtor(k,2,i)*ddnewtor(l,2,j)) + v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) & + +ccnewtor(k,1,i)*ddnewtor(l,2,j)) +#else + v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) & + -ccnewtor(k,2,i)*ddnewtor(l,2,j)) + v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) & + +ccnewtor(k,1,i)*ddnewtor(l,2,j)) +#endif + enddo + enddo +!c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma) + enddo + enddo +#else + write (iout,*) "TOR_MODE>1 only with NEWCORR" + stop +#endif + ENDIF ! TOR_MODE + + if (tor_mode.gt.0 .and. lprint) then +!c Print valence-torsional parameters + write (iout,'(a)') & + "Parameters of the valence-torsional potentials" + do i=-ntortyp+1,ntortyp-1 + do j=-ntortyp+1,ntortyp-1 + write (iout,'(3a)') "Type ",toronelet(i),toronelet(j) + write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc" + do k=1,nterm_kcc(j,i) + do l=1,nterm_kcc_Tb(j,i) + do ll=1,nterm_kcc_Tb(j,i) + write (iout,'(3i5,2f15.4)')& + k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i) + enddo + enddo + enddo + enddo + enddo + endif #endif allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1) read (itorp_nucl,*,end=113,err=113) ntortyp_nucl @@ -1930,139 +2500,6 @@ ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local ! interaction energy of the Gly, Ala, and Pro prototypes. ! - - if (lprint) then - write (iout,*) - write (iout,*) "Coefficients of the cumulants" - endif - read (ifourier,*) nloctyp -!write(iout,*) "nloctyp",nloctyp -!el from module energy------- - allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor) - allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor) - allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor) - allocate(cc(2,2,-nloctyp-1:nloctyp+1)) - allocate(dd(2,2,-nloctyp-1:nloctyp+1)) - allocate(ee(2,2,-nloctyp-1:nloctyp+1)) - allocate(ctilde(2,2,-nloctyp-1:nloctyp+1)) - allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor) -! el - b1(1,:)=0.0d0 - b1(2,:)=0.0d0 -!-------------------------------- - - do i=0,nloctyp-1 - read (ifourier,*,end=115,err=115) - read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13) - if (lprint) then - write (iout,*) 'Type',i - write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13) - endif - B1(1,i) = buse(3) - B1(2,i) = buse(5) - B1(1,-i) = buse(3) - B1(2,-i) = -buse(5) -! buse1(1,i)=0.0d0 -! buse1(2,i)=0.0d0 - B1tilde(1,i) = buse(3) - B1tilde(2,i) =-buse(5) - B1tilde(1,-i) =-buse(3) - B1tilde(2,-i) =buse(5) -! buse1tilde(1,i)=0.0d0 -! buse1tilde(2,i)=0.0d0 - B2(1,i) = buse(2) - B2(2,i) = buse(4) - B2(1,-i) =buse(2) - B2(2,-i) =-buse(4) - -! buse2(1,i)=0.0d0 -! buse2(2,i)=0.0d0 - CC(1,1,i)= buse(7) - CC(2,2,i)=-buse(7) - CC(2,1,i)= buse(9) - CC(1,2,i)= buse(9) - CC(1,1,-i)= buse(7) - CC(2,2,-i)=-buse(7) - CC(2,1,-i)=-buse(9) - CC(1,2,-i)=-buse(9) -! CC(1,1,i)=0.0d0 -! CC(2,2,i)=0.0d0 -! CC(2,1,i)=0.0d0 -! CC(1,2,i)=0.0d0 - Ctilde(1,1,i)=buse(7) - Ctilde(1,2,i)=buse(9) - Ctilde(2,1,i)=-buse(9) - Ctilde(2,2,i)=buse(7) - Ctilde(1,1,-i)=buse(7) - Ctilde(1,2,-i)=-buse(9) - Ctilde(2,1,-i)=buse(9) - Ctilde(2,2,-i)=buse(7) - -! Ctilde(1,1,i)=0.0d0 -! Ctilde(1,2,i)=0.0d0 -! Ctilde(2,1,i)=0.0d0 -! Ctilde(2,2,i)=0.0d0 - DD(1,1,i)= buse(6) - DD(2,2,i)=-buse(6) - DD(2,1,i)= buse(8) - DD(1,2,i)= buse(8) - DD(1,1,-i)= buse(6) - DD(2,2,-i)=-buse(6) - DD(2,1,-i)=-buse(8) - DD(1,2,-i)=-buse(8) -! DD(1,1,i)=0.0d0 -! DD(2,2,i)=0.0d0 -! DD(2,1,i)=0.0d0 -! DD(1,2,i)=0.0d0 - Dtilde(1,1,i)=buse(6) - Dtilde(1,2,i)=buse(8) - Dtilde(2,1,i)=-buse(8) - Dtilde(2,2,i)=buse(6) - Dtilde(1,1,-i)=buse(6) - Dtilde(1,2,-i)=-buse(8) - Dtilde(2,1,-i)=buse(8) - Dtilde(2,2,-i)=buse(6) - -! Dtilde(1,1,i)=0.0d0 -! Dtilde(1,2,i)=0.0d0 -! Dtilde(2,1,i)=0.0d0 -! Dtilde(2,2,i)=0.0d0 - EE(1,1,i)= buse(10)+buse(11) - EE(2,2,i)=-buse(10)+buse(11) - EE(2,1,i)= buse(12)-buse(13) - EE(1,2,i)= buse(12)+buse(13) - EE(1,1,-i)= buse(10)+buse(11) - EE(2,2,-i)=-buse(10)+buse(11) - EE(2,1,-i)=-buse(12)+buse(13) - EE(1,2,-i)=-buse(12)-buse(13) - -! ee(1,1,i)=1.0d0 -! ee(2,2,i)=1.0d0 -! ee(2,1,i)=0.0d0 -! ee(1,2,i)=0.0d0 -! ee(2,1,i)=ee(1,2,i) - enddo - if (lprint) then - do i=1,nloctyp - write (iout,*) 'Type',i - write (iout,*) 'B1' - write(iout,*) B1(1,i),B1(2,i) - write (iout,*) 'B2' - write(iout,*) B2(1,i),B2(2,i) - write (iout,*) 'CC' - do j=1,2 - write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i) - enddo - write(iout,*) 'DD' - do j=1,2 - write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i) - enddo - write(iout,*) 'EE' - do j=1,2 - write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i) - enddo - enddo - endif ! ! Read electrostatic-interaction parameters ! @@ -2815,6 +3252,8 @@ 118 write (iout,*) "Error reading SCp interaction parameters." goto 999 119 write (iout,*) "Error reading SCCOR parameters" + go to 999 + 121 write (iout,*) "Error in Czybyshev parameters" 999 continue #ifdef MPI call MPI_Finalize(Ierror) @@ -3759,6 +4198,10 @@ call readi(controlcard,'SHIELD',shield_mode,0) !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then write(iout,*) "shield_mode",shield_mode + call readi(controlcard,'TORMODE',tor_mode,0) +!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) "torsional and valence angle mode",tor_mode + !C Varibles set size of box with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 protein=index(controlcard,"PROTEIN").gt.0 diff --git a/source/unres/unres.F90 b/source/unres/unres.F90 index 3713a37..bc9da96 100644 --- a/source/unres/unres.F90 +++ b/source/unres/unres.F90 @@ -77,6 +77,8 @@ ! Read force field parameters and job setup data call readrtns call flush(iout) + write (iout,*) "After readrtns" + call cartprint ! if (me.eq.king .or. .not. out1file) then write (iout,'(2a/)') &