diff --git a/examples/amoeba/data.ubiquitin b/examples/amoeba/data.ubiquitin index a7b8684618..05a0673ce8 100644 --- a/examples/amoeba/data.ubiquitin +++ b/examples/amoeba/data.ubiquitin @@ -13,6 +13,8 @@ LAMMPS data file created from Tinker ubiquitin.xyz and amoeba_ubiquitin.prm file 0 54.99 xlo xhi 0 41.91 ylo yhi 0 41.91 zlo zhi +6 pitorsion types +106 pitorsions Masses @@ -26346,6 +26348,124 @@ UreyBradley Coeffs 110 0.0 0.0 111 -7.6 1.5537 +PiTorsions + +1 1 2 4 3 20 21 24 +2 1 21 23 22 37 38 41 +3 1 27 29 28 30 35 36 +4 1 38 40 39 56 57 60 +5 1 57 59 58 76 77 80 +6 2 62 65 63 64 66 71 +7 2 62 64 63 65 67 72 +8 2 63 71 64 66 68 73 +9 2 63 72 65 67 68 74 +10 2 64 73 66 68 67 75 +11 2 65 74 67 68 66 75 +12 1 77 79 78 92 93 96 +13 1 93 95 94 114 115 118 +14 1 115 117 116 128 129 132 +15 1 129 131 130 147 148 151 +16 1 148 150 149 161 162 165 +17 1 162 164 163 168 169 172 +18 1 169 171 170 190 191 194 +19 1 191 193 192 204 205 208 +20 1 205 207 206 223 224 227 +21 1 224 226 225 237 238 241 +22 1 238 240 239 256 257 260 +23 1 257 259 258 271 272 275 +24 1 272 274 273 287 288 291 +25 3 288 290 289 302 303 309 +26 1 303 305 304 316 317 320 +27 1 317 319 318 327 328 331 +28 1 328 330 329 339 340 343 +29 1 340 342 341 353 354 357 +30 1 354 356 355 372 373 376 +31 1 373 375 374 387 388 391 +32 1 388 390 389 401 402 405 +33 1 393 395 394 396 399 400 +34 1 402 404 403 417 418 421 +35 1 418 420 419 439 440 443 +36 1 440 442 441 449 450 453 +37 1 450 452 451 471 472 475 +38 1 472 474 473 490 491 494 +39 1 491 493 492 507 508 511 +40 1 497 499 498 500 505 506 +41 1 508 510 509 519 520 523 +42 1 520 522 521 541 542 545 +43 1 542 544 543 556 557 560 +44 1 557 559 558 563 564 567 +45 3 564 566 565 582 583 589 +46 3 583 585 584 596 597 603 +47 1 597 599 598 610 611 614 +48 1 611 613 612 622 623 626 +49 1 623 625 624 639 640 643 +50 1 629 631 630 632 637 638 +51 1 640 642 641 656 657 660 +52 1 646 648 647 649 654 655 +53 1 657 659 658 680 681 684 +54 1 681 683 682 699 700 703 +55 1 700 702 701 718 719 722 +56 1 719 721 720 738 739 742 +57 2 724 727 725 726 728 733 +58 2 724 726 725 727 729 734 +59 2 725 733 726 728 730 735 +60 2 725 734 727 729 730 736 +61 2 726 735 728 730 729 737 +62 2 727 736 729 730 728 737 +63 1 739 741 740 748 749 752 +64 1 749 751 750 755 756 759 +65 1 756 758 757 777 778 781 +66 1 778 780 779 794 795 798 +67 1 784 786 785 787 792 793 +68 1 795 797 796 813 814 817 +69 1 814 816 815 828 829 832 +70 1 829 831 830 840 841 844 +71 1 841 843 842 847 848 851 +72 1 848 850 849 871 872 875 +73 1 872 874 873 885 886 889 +74 1 886 888 887 904 905 908 +75 1 905 907 906 915 916 919 +76 1 916 918 917 927 928 931 +77 1 928 930 929 948 949 952 +78 2 933 936 934 935 937 943 +79 2 933 935 934 936 938 944 +80 2 934 943 935 937 939 945 +81 2 934 944 936 938 939 946 +82 2 935 945 937 939 938 940 +83 2 936 946 938 939 937 940 +84 1 949 951 950 962 963 966 +85 1 954 956 955 957 960 961 +86 1 963 965 964 981 982 985 +87 1 982 984 983 998 999 1002 +88 1 988 990 989 991 996 997 +89 1 999 1001 1000 1020 1021 1024 +90 1 1021 1023 1022 1035 1036 1039 +91 1 1036 1038 1037 1046 1047 1050 +92 1 1047 1049 1048 1060 1061 1064 +93 1 1061 1063 1062 1079 1080 1083 +94 1 1080 1082 1081 1097 1098 1101 +95 4 1085 1088 1086 1087 1089 1093 +96 5 1085 1087 1086 1088 1090 1094 +97 6 1086 1093 1087 1089 1090 1095 +98 4 1086 1094 1088 1090 1089 1096 +99 6 1087 1095 1089 1090 1088 1096 +100 1 1098 1100 1099 1116 1117 1120 +101 1 1117 1119 1118 1132 1133 1136 +102 1 1133 1135 1134 1151 1152 1155 +103 1 1152 1154 1153 1175 1176 1179 +104 1 1176 1178 1177 1194 1195 1198 +105 1 1195 1197 1196 1218 1219 1222 +106 1 1219 1221 1220 1225 1226 1229 + +PiTorsion Coeffs + +1 6.85 +2 6.85 +3 6.85 +4 6.85 +5 6.85 +6 6.85 + Tinker Types 1 234 diff --git a/examples/amoeba/in.ubiquitin b/examples/amoeba/in.ubiquitin index 0d7a708826..367030ba74 100644 --- a/examples/amoeba/in.ubiquitin +++ b/examples/amoeba/in.ubiquitin @@ -14,15 +14,16 @@ improper_style amoeba # per-atom properties required by AMOEBA or HIPPO fix amtype all property/atom i_amtype ghost yes -#fix pitorsion all pitorsion +fix pitorsion all pitorsion fix extra all property/atom & i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes fix extra2 all property/atom i_polaxe -read_data data.ubiquitin fix amtype NULL "Tinker Types" -#read_data data.ubiquitin fix amtype NULL "Tinker Types" & -# fix pitorsion pitorsions PiTorsions +#read_data data.ubiquitin fix amtype NULL "Tinker Types" +read_data data.ubiquitin fix amtype NULL "Tinker Types" & + fix pitorsion "pitorsion types" "PiTorsion Coeffs" & + fix pitorsion pitorsions PiTorsions pair_style amoeba pair_coeff * * amoeba_ubiquitin.prm amoeba_ubiquitin.key diff --git a/src/AMOEBA/angle_amoeba.cpp b/src/AMOEBA/angle_amoeba.cpp index a578396d10..e81bbff00f 100644 --- a/src/AMOEBA/angle_amoeba.cpp +++ b/src/AMOEBA/angle_amoeba.cpp @@ -103,10 +103,6 @@ void AngleAmoeba::compute(int eflag, int vflag) ev_init(eflag,vflag); - // DEBUG - int ubn = 0; - double eub = 0.0; - for (n = 0; n < nanglelist; n++) { i1 = anglelist[n][0]; i2 = anglelist[n][1]; @@ -131,15 +127,8 @@ void AngleAmoeba::compute(int eflag, int vflag) // Urey-Bradley H-H bond term within water molecules uflag = ubflag[type]; - - if (uflag) { - ubn++; - eub += tinker_urey_bradley(i1,i3,type,eflag); - } + if (uflag) tinker_urey_bradley(i1,i3,type,eflag); } - - // DEBUG - printf("UreyBradley n %d eng %g\n",ubn,eub); } /* ---------------------------------------------------------------------- */ @@ -543,7 +532,7 @@ void AngleAmoeba::tinker_bondangle(int i1, int i2, int i3, int type, int eflag) /* ---------------------------------------------------------------------- */ -double AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag) +void AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag) { double delx,dely,delz; double rsq,r,dr,rk; @@ -585,8 +574,6 @@ double AngleAmoeba::tinker_urey_bradley(int i1, int i2, int type, int eflag) } if (evflag) ev_tally2(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz); - - return ebond; } /* ---------------------------------------------------------------------- */ diff --git a/src/AMOEBA/angle_amoeba.h b/src/AMOEBA/angle_amoeba.h index 281f4a9adb..69a782b48b 100644 --- a/src/AMOEBA/angle_amoeba.h +++ b/src/AMOEBA/angle_amoeba.h @@ -46,7 +46,7 @@ class AngleAmoeba : public Angle { void tinker_angle(int, int, int, int, int); void tinker_anglep(int, int, int, int, int); void tinker_bondangle(int, int, int, int, int); - double tinker_urey_bradley(int, int, int, int); + void tinker_urey_bradley(int, int, int, int); void allocate(); }; diff --git a/src/AMOEBA/fix_pitorsion.cpp b/src/AMOEBA/fix_pitorsion.cpp index ec8f0d71be..939220d95d 100644 --- a/src/AMOEBA/fix_pitorsion.cpp +++ b/src/AMOEBA/fix_pitorsion.cpp @@ -89,6 +89,10 @@ FixPiTorsion::FixPiTorsion(LAMMPS *lmp, int narg, char **arg) : npitorsion_list = 0; max_pitorsion_list = 0; pitorsion_list = nullptr; + + // pitorsion coeff + + kpit = nullptr; } /* --------------------------------------------------------------------- */ @@ -114,6 +118,10 @@ FixPiTorsion::~FixPiTorsion() // compute list memory->destroy(pitorsion_list); + + // coeff + + memory->destroy(kpit); } /* ---------------------------------------------------------------------- */ @@ -185,19 +193,19 @@ void FixPiTorsion::pre_neighbor() int i,m,atom1,atom2,atom3,atom4,atom5,atom6; // guesstimate initial length of local pitorsion list - // if npitorsion_global was not set (due to read_restart, no read_data), + // if npitorsions was not set (due to read_restart, no read_data), // then list will grow by LISTDELTA chunks if (max_pitorsion_list == 0) { - if (nprocs == 1) max_pitorsion_list = npitorsion_global; + if (nprocs == 1) max_pitorsion_list = npitorsions; else max_pitorsion_list = - static_cast (LB_FACTOR*npitorsion_global/nprocs); + static_cast (LB_FACTOR*npitorsions/nprocs); memory->create(pitorsion_list,max_pitorsion_list,7, "pitorsion:pitorsion_list"); } int nlocal = atom->nlocal; - pitorsion_list = 0; + npitorsion_list = 0; for (i = 0; i < nlocal; i++) { for (m = 0; m < num_pitorsion[i]; m++) { @@ -209,7 +217,7 @@ void FixPiTorsion::pre_neighbor() atom6 = atom->map(pitorsion_atom6[i][m]); if (atom1 == -1 || atom2 == -1 || atom3 == -1 || - atom4 == -1 || atom5 == -1) + atom4 == -1 || atom5 == -1 || atom6 == -1) error->one(FLERR,"PiTorsion atoms {} {} {} {} {} {} missing on " "proc {} at step {}", pitorsion_atom1[i][m],pitorsion_atom2[i][m], @@ -531,7 +539,9 @@ double FixPiTorsion::compute_scalar() void FixPiTorsion::read_data_header(char *line) { if (strstr(line,"pitorsions")) { - sscanf(line,BIGINT_FORMAT,&npitorsion_global); + sscanf(line,BIGINT_FORMAT,&npitorsions); + } else if (strstr(line,"pitorsion types")) { + sscanf(line,"%d",&npitorsion_types); } else error->all(FLERR,"Invalid read data header line for fix pitorsion"); // didn't set in constructor because this fix could be defined @@ -549,123 +559,172 @@ void FixPiTorsion::read_data_header(char *line) void FixPiTorsion::read_data_section(char *keyword, int n, char *buf, tagint id_offset) { - int m,tmp,itype; - tagint atom1,atom2,atom3,atom4,atom5,atom6; - char *next; + int which; - next = strchr(buf,'\n'); - *next = '\0'; - int nwords = utils::count_words(utils::trim_comment(buf)); - *next = '\n'; + if (strstr(keyword,"PiTorsions")) { + sscanf(keyword,BIGINT_FORMAT,&npitorsions); + which = 0; + } else if (strstr(keyword,"PiTorsion Coeffs")) { + sscanf(keyword,"%d",&npitorsion_types); + which = 1; + } else error->all(FLERR,"Invalid read data section for fix pitorsion"); + + // loop over lines of PiTorsion Coeffs + // tokenize the line into values + // initialize kpit vector + + if (which == 1) { + int itype; + double value; + char *next; - if (nwords != 7) - error->all(FLERR,"Incorrect {} format in data file",keyword); + memory->create(kpit,npitorsion_types+1,"pitorsion:kpit"); + + for (int i = 0; i < n; i++) { + next = strchr(buf,'\n'); + *next = '\0'; + sscanf(buf,"%d %g",&itype,&value); + if (itype <= 0 || itype > npitorsion_types) + error->all(FLERR,"Incorrect args for pitorsion coefficients"); + kpit[itype] = value; + buf = next + 1; + } + } // loop over lines of PiTorsions // tokenize the line into values // add PiTorsion to one or more of my atoms, depending on newton_bond - for (int i = 0; i < n; i++) { + if (which == 0) { + + int m,tmp,itype; + tagint atom1,atom2,atom3,atom4,atom5,atom6; + char *next; + next = strchr(buf,'\n'); *next = '\0'; - sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT - " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, - &tmp,&itype,&atom1,&atom2,&atom3,&atom4,&atom5,&atom6); + int nwords = utils::count_words(utils::trim_comment(buf)); + *next = '\n'; + + if (nwords != 8) + error->all(FLERR,"Incorrect {} format in data file",keyword); - atom1 += id_offset; - atom2 += id_offset; - atom3 += id_offset; - atom4 += id_offset; - atom5 += id_offset; - atom6 += id_offset; + for (int i = 0; i < n; i++) { + next = strchr(buf,'\n'); + *next = '\0'; + sscanf(buf,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT + " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT, + &tmp,&itype,&atom1,&atom2,&atom3,&atom4,&atom5,&atom6); - if ((m = atom->map(atom1)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; + atom1 += id_offset; + atom2 += id_offset; + atom3 += id_offset; + atom4 += id_offset; + atom5 += id_offset; + atom6 += id_offset; + + // atom3 owns the pitorsion when newton_bond on + + if ((m = atom->map(atom3)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; + } + + if (newton_bond == 0) { + if ((m = atom->map(atom1)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; + } + } + + if (newton_bond == 0) { + if ((m = atom->map(atom2)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; + } + } + + if (newton_bond == 0) { + if ((m = atom->map(atom4)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; + } + } + + if (newton_bond == 0) { + if ((m = atom->map(atom5)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; + } + } + + if (newton_bond == 0) { + if ((m = atom->map(atom6)) >= 0) { + if (num_pitorsion[m] == PITORSIONMAX) + error->one(FLERR,"Too many PiTorsions for one atom"); + pitorsion_type[m][num_pitorsion[m]] = itype; + pitorsion_atom1[m][num_pitorsion[m]] = atom1; + pitorsion_atom2[m][num_pitorsion[m]] = atom2; + pitorsion_atom3[m][num_pitorsion[m]] = atom3; + pitorsion_atom4[m][num_pitorsion[m]] = atom4; + pitorsion_atom5[m][num_pitorsion[m]] = atom5; + pitorsion_atom6[m][num_pitorsion[m]] = atom6; + num_pitorsion[m]++; + } + } + + buf = next + 1; } - - if ((m = atom->map(atom2)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } - - if ((m = atom->map(atom3)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } - - if ((m = atom->map(atom4)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } - - if ((m = atom->map(atom5)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } - - if ((m = atom->map(atom6)) >= 0) { - if (num_pitorsion[m] == PITORSIONMAX) - error->one(FLERR,"Too many PiTorsions for one atom"); - pitorsion_type[m][num_pitorsion[m]] = itype; - pitorsion_atom1[m][num_pitorsion[m]] = atom1; - pitorsion_atom2[m][num_pitorsion[m]] = atom2; - pitorsion_atom3[m][num_pitorsion[m]] = atom3; - pitorsion_atom4[m][num_pitorsion[m]] = atom4; - pitorsion_atom5[m][num_pitorsion[m]] = atom5; - pitorsion_atom6[m][num_pitorsion[m]] = atom6; - num_pitorsion[m]++; - } - - buf = next + 1; } } /* ---------------------------------------------------------------------- */ -bigint FixPiTorsion::read_data_skip_lines(char * /*keyword*/) +bigint FixPiTorsion::read_data_skip_lines(char *keyword) { - return npitorsion_global; + if (strcmp(keyword,"PiTorsions") == 0) return npitorsions; + if (strcmp(keyword,"PiTorsion Coeffs") == 0) return (bigint) npitorsion_types; + return 0; } /* ---------------------------------------------------------------------- @@ -673,9 +732,11 @@ bigint FixPiTorsion::read_data_skip_lines(char * /*keyword*/) only called by proc 0 ------------------------------------------------------------------------- */ -void FixPiTorsion::write_data_header(FILE *fp, int /*mth*/) +void FixPiTorsion::write_data_header(FILE *fp, int mth) { - fprintf(fp,BIGINT_FORMAT " pitorsions\n",npitorsion_global); + if (mth == 0) fprintf(fp,BIGINT_FORMAT " pitorsions\n",npitorsions); + else if (mth == 1) + fprintf(fp,BIGINT_FORMAT " pitorsion types\n",npitorsion_types); } /* ---------------------------------------------------------------------- @@ -686,7 +747,7 @@ void FixPiTorsion::write_data_header(FILE *fp, int /*mth*/) ny = columns = type + 6 atom IDs ------------------------------------------------------------------------- */ -void FixPiTorsion::write_data_section_size(int /*mth*/, int &nx, int &ny) +void FixPiTorsion::write_data_section_size(int mth, int &nx, int &ny) { int i,m; @@ -708,7 +769,7 @@ void FixPiTorsion::write_data_section_size(int /*mth*/, int &nx, int &ny) buf allocated by caller as owned PiTorsions by 7 ------------------------------------------------------------------------- */ -void FixPiTorsion::write_data_section_pack(int /*mth*/, double **buf) +void FixPiTorsion::write_data_section_pack(int mth, double **buf) { int i,m; @@ -741,9 +802,10 @@ void FixPiTorsion::write_data_section_pack(int /*mth*/, double **buf) only called by proc 0 ------------------------------------------------------------------------- */ -void FixPiTorsion::write_data_section_keyword(int /*mth*/, FILE *fp) +void FixPiTorsion::write_data_section_keyword(int mth, FILE *fp) { - fprintf(fp,"\nPiTorsion\n\n"); + if (mth == 0) fprintf(fp,"\nPiTorsions\n\n"); + else if (mth == 1) fprintf(fp,"\nPiTorsion Coeffs\n\n"); } /* ---------------------------------------------------------------------- @@ -753,7 +815,7 @@ void FixPiTorsion::write_data_section_keyword(int /*mth*/, FILE *fp) only called by proc 0 ------------------------------------------------------------------------- */ -void FixPiTorsion::write_data_section(int /*mth*/, FILE *fp, +void FixPiTorsion::write_data_section(int mth, FILE *fp, int n, double **buf, int index) { for (int i = 0; i < n; i++) @@ -779,7 +841,7 @@ void FixPiTorsion::write_restart(FILE *fp) if (comm->me == 0) { int size = sizeof(bigint); fwrite(&size,sizeof(int),1,fp); - fwrite(&npitorsion_global,sizeof(bigint),1,fp); + fwrite(&npitorsions,sizeof(bigint),1,fp); } } @@ -789,7 +851,7 @@ void FixPiTorsion::write_restart(FILE *fp) void FixPiTorsion::restart(char *buf) { - npitorsion_global = *((bigint *) buf); + npitorsions = *((bigint *) buf); } /* ---------------------------------------------------------------------- diff --git a/src/AMOEBA/fix_pitorsion.h b/src/AMOEBA/fix_pitorsion.h index efb448a408..7e3258c66a 100644 --- a/src/AMOEBA/fix_pitorsion.h +++ b/src/AMOEBA/fix_pitorsion.h @@ -68,7 +68,8 @@ class FixPiTorsion : public Fix { int nprocs, me; int newton_bond, eflag_caller; int ilevel_respa; - bigint npitorsion_global; + bigint npitorsions; + int npitorsion_types; double epitorsion; double *kpit; diff --git a/src/modify.cpp b/src/modify.cpp index eff6ea8337..ef738ecc67 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -810,7 +810,7 @@ Fix *Modify::add_fix(int narg, char **arg, int trysuffix) const char *exceptions[] = {"GPU", "OMP", "INTEL", "property/atom", "cmap", "cmap3", "rx", - "deprecated", "STORE/KIM", nullptr}; + "deprecated", "STORE/KIM", "pitorsion", nullptr}; if (domain->box_exist == 0) { int m; diff --git a/src/read_data.cpp b/src/read_data.cpp index 5b68f9d94e..25e2db39bf 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -729,7 +729,6 @@ void ReadData::command(int narg, char **arg) if (firstpass) fix(fix_index[i],keyword); else skip_lines(modify->fix[fix_index[i]]-> read_data_skip_lines(keyword)); - parse_keyword(0); break; } if (i == nfix) diff --git a/tools/tinker/tinker2lmp.py b/tools/tinker/tinker2lmp.py index 95e0013e91..dfacebd832 100644 --- a/tools/tinker/tinker2lmp.py +++ b/tools/tinker/tinker2lmp.py @@ -718,7 +718,7 @@ for atom1 in id: # DEBUG - if uncommented, no PiTorsions appear in data file -pitorsionlist = [] +#pitorsionlist = [] # ---------------------------------------- # create lists of bond/angle/dihedral/improper types @@ -1055,32 +1055,31 @@ for m,params in enumerate(prm.pitorsionparams): flags = len(prm.pitorsionparams)*[0] pitorsiontype = [] pitorsionparams = [] +pitorsionlist_reduced = [] for tmp1,tmp2,atom1,atom2,tmp3,tmp4 in pitorsionlist: type1 = type[atom1-1] type2 = type[atom2-1] c1 = classes[type1-1] c2 = classes[type2-1] - - if (c1,c2) in pitdict: m,params = pitdict[(c1,c2)] - elif (c2,c1) in pitdict: m,params = pitdict[(c2,c1)] - else: - # NOTE: probably this PiT should be deleted from enumeration list - # similar to olist for impropers above - # IMPORTANT: look at other errors and messages t2lmp.py is producing - error("PiTorsion not found: %d %d: %d %d" % (atom1,atom2,c1,c2)) - pitorsiontype.append(0) - continue - - if not flags[m]: - v1 = params[2:] - pitorsionparams.append(v1) - flags[m] = len(pitorsionparams) - pitorsiontype.append(flags[m]) -print "PRM pitor param",len(prm.pitorsionparams) -print "PiTor list",len(pitorsionlist) -#print "Flags",flags + # 6-tuple is only a PiTorsion if central 2 atoms math an entry in PRM file + # pitorsionlist_reduced = list of just these 6-tuples + + if (c1,c2) in pitdict or (c2,c1) in pitdict: + if (c1,c2) in pitdict: m,params = pitdict[(c1,c2)] + else: m,params = pitdict[(c2,c1)] + pitorsionlist_reduced.append((tmp1,tmp2,atom1,atom2,tmp3,tmp4)) + + if not flags[m]: + v1 = params[2:] + pitorsionparams.append(v1) + flags[m] = len(pitorsionparams) + pitorsiontype.append(flags[m]) + +# replace original pitorsionlist with reduced version + +pitorsionlist = pitorsionlist_reduced # ---------------------------------------- # assign each atom to a Tinker group