ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

initial version of fix pitorsion

Browse Source
This commit is contained in:
Steve Plimpton
2022-03-09 15:37:19 -07:00
parent 844ea0ab8e
commit bbe065e649
9 changed files with 326 additions and 157 deletions
Download Patch File Download Diff File Expand all files Collapse all files

120
examples/amoeba/data.ubiquitin
View File

@ -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

9
examples/amoeba/in.ubiquitin
View File

@ -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

17
src/AMOEBA/angle_amoeba.cpp
View File

@ -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,17 +127,10 @@ 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);
}
/* ---------------------------------------------------------------------- */
void AngleAmoeba::tinker_angle(int i1, int i2, int i3, int type, int eflag)
@ -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;
}
/* ---------------------------------------------------------------------- */

2
src/AMOEBA/angle_amoeba.h
View File

@ -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();
};

156
src/AMOEBA/fix_pitorsion.cpp
View File

@ -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<int> (LB_FACTOR*npitorsion_global/nprocs);
static_cast<int> (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,6 +559,44 @@ void FixPiTorsion::read_data_header(char *line)
void FixPiTorsion::read_data_section(char *keyword, int n, char *buf,
tagint id_offset)
{
int which;
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;
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
if (which == 0) {
int m,tmp,itype;
tagint atom1,atom2,atom3,atom4,atom5,atom6;
char *next;
@ -558,13 +606,9 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf,
int nwords = utils::count_words(utils::trim_comment(buf));
*next = '\n';
if (nwords != 7)
if (nwords != 8)
error->all(FLERR,"Incorrect {} format in data file",keyword);
// 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++) {
next = strchr(buf,'\n');
*next = '\0';
@ -579,31 +623,7 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf,
atom5 += id_offset;
atom6 += id_offset;
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 ((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]++;
}
// atom3 owns the pitorsion when newton_bond on
if ((m = atom->map(atom3)) >= 0) {
if (num_pitorsion[m] == PITORSIONMAX)
@ -618,6 +638,37 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf,
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");
@ -630,7 +681,9 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf,
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");
@ -643,7 +696,9 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf,
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");
@ -656,16 +711,20 @@ void FixPiTorsion::read_data_section(char *keyword, int n, char *buf,
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);
}
/* ----------------------------------------------------------------------

3
src/AMOEBA/fix_pitorsion.h
View File

@ -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;

2
src/modify.cpp
View File

@ -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;

1
src/read_data.cpp
View File

@ -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)

23
tools/tinker/tinker2lmp.py
View File

@ -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,6 +1055,7 @@ 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]
@ -1062,15 +1063,13 @@ for tmp1,tmp2,atom1,atom2,tmp3,tmp4 in pitorsionlist:
c1 = classes[type1-1]
c2 = classes[type2-1]
# 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)]
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
else: m,params = pitdict[(c2,c1)]
pitorsionlist_reduced.append((tmp1,tmp2,atom1,atom2,tmp3,tmp4))
if not flags[m]:
v1 = params[2:]
@ -1078,9 +1077,9 @@ for tmp1,tmp2,atom1,atom2,tmp3,tmp4 in pitorsionlist:
flags[m] = len(pitorsionparams)
pitorsiontype.append(flags[m])
print "PRM pitor param",len(prm.pitorsionparams)
print "PiTor list",len(pitorsionlist)
#print "Flags",flags
# replace original pitorsionlist with reduced version
pitorsionlist = pitorsionlist_reduced
# ----------------------------------------
# assign each atom to a Tinker group