From: Berk Hess Date: Sun, 15 Sep 2019 21:01:43 +0000 (+0200) Subject: Add virtual site type 2FD X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=b766e6ebfa539173b4b3580ca5356533a27a8caa;p=alexxy%2Fgromacs.git Add virtual site type 2FD Fixes #2451 Change-Id: Ide9dcd829284567010435ed5cfb55194aed12dcf --- diff --git a/docs/reference-manual/functions/interaction-methods.rst b/docs/reference-manual/functions/interaction-methods.rst index 45812f98bc..c036c671df 100644 --- a/docs/reference-manual/functions/interaction-methods.rst +++ b/docs/reference-manual/functions/interaction-methods.rst @@ -170,7 +170,7 @@ However, in the general case redistribution should be done first. .. figure:: plots/dummies.* :width: 15.00000cm - The six different types of virtual site construction in . The + The seven different types of virtual site construction. The constructing atoms are shown as black circles, the virtual sites in gray. @@ -206,6 +206,30 @@ can be constructed from “particles” that are simpler virtual sites. - In this case the virtual site is on the line through atoms :math:`i` and :math:`j`. +- On the line through two atoms, with a fixed distance + (:numref:`Fig. %s ` 2fd): + + .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + a \frac{ \mathbf{r}_ij } + { | \mathbf{r}_ij | } + :label: eqnvsite2fdatom + +- In this case the virtual site is on the line through the the other two + particles at a distance of :math:`|a|` from :math:`i`. The force on + particles :math:`i` and :math:`j` due to the force on the virtual site + can be computed as: + + .. math:: \begin{array}{lcr} + \mathbf{F}_i &=& \displaystyle \mathbf{F}_{s} - \gamma ( \mathbf{F}_is - \mathbf{p} ) \\[1ex] + \mathbf{F}_j &=& \displaystyle \gamma (\mathbf{F}_{s} - \mathbf{p}) \\[1ex] + \end{array} + ~\mbox{~ where~ }~ + \begin{array}{c} + \displaystyle \gamma = \frac{a}{ | \mathbf{r}_ij | } \\[2ex] + \displaystyle \mathbf{p} = \frac{ \mathbf{r}_{is} \cdot \mathbf{F}_{s} } + { \mathbf{r}_{is} \cdot \mathbf{r}_is } \mathbf{r}_is + \end{array} + :label: eqnvsite2fdforce + - As a linear combination of three atoms (:numref:`Fig. %s ` 3): @@ -218,8 +242,8 @@ can be constructed from “particles” that are simpler virtual sites. - In the plane of three atoms, with a fixed distance (:numref:`Fig. %s ` 3fd): - .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + b \frac{ \mathbf{r}_ij + a \mathbf{r}_{jk} } - { | \mathbf{r}_ij + a \mathbf{r}_{jk} | } + .. math:: \mathbf{r}_s ~=~ \mathbf{r}_i + b \frac{ (1 - a) \mathbf{r}_ij + a \mathbf{r}_{jk} } + { | (1 - a) \mathbf{r}_ij + a \mathbf{r}_{jk} | } :label: eqnvsiteplane3atom - In this case the virtual site is in the plane of the other three diff --git a/docs/reference-manual/functions/plots/dummies.fig b/docs/reference-manual/functions/plots/dummies.fig index f868c8a431..074884422c 100644 --- a/docs/reference-manual/functions/plots/dummies.fig +++ b/docs/reference-manual/functions/plots/dummies.fig @@ -1,105 +1,119 @@ -#FIG 3.1 +#FIG 3.2 Produced by xfig version 3.2.7a Portrait Center Inches +A4 +100.00 +Single +-2 1200 2 5 1 2 1 -1 -1 0 0 -1 3.000 0 1 0 0 11851.786 1198.214 12525 1200 11850 525 11550 1800 -1 3 0 1 -1 7 0 0 0 0.000 1 0.0000 6825 1200 150 150 6825 1200 6975 1350 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 6825 3600 150 150 6825 3600 6975 3600 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 7500 2400 150 150 7500 2400 7650 2400 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 11175 2400 150 150 11175 2400 11325 2400 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 9825 2400 150 150 9825 2400 9975 2400 -1 3 2 1 -1 -1 0 0 -1 4.000 1 0.0000 6825 2400 75 75 6825 2400 6900 2475 -1 3 0 1 -1 7 0 0 0 0.000 1 0.0000 4050 1200 150 150 4050 1200 4200 1350 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 4050 3600 150 150 4050 3600 4200 3600 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 4725 2400 150 150 4725 2400 4875 2400 1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 2925 2400 150 150 2925 2400 3075 2400 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 1125 2400 150 150 1125 2400 1275 2550 1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 2025 2400 150 150 2025 2400 2175 2400 -1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 5625 2400 168 168 5625 2400 5775 2475 -1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 8400 2400 150 150 8400 2400 8550 2400 -1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 12750 1200 150 150 12750 1200 12900 1200 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 11850 1200 150 150 11850 1200 12000 1200 -1 3 0 1 -1 -1 0 0 20 0.000 1 0.0000 13740 2355 150 150 13740 2355 13890 2355 -1 3 0 1 7 -1 2 0 20 0.000 1 0.0000 15165 2730 150 150 15165 2730 15315 2880 -1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 15390 3180 168 168 15390 3180 15540 3255 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 15165 1980 150 150 15165 1980 15315 2130 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 16800 1593 168 168 16800 1593 16632 1425 -1 3 0 1 -1 -1 2 0 62 0.000 1 0.0000 17100 2400 150 150 17100 2400 17250 2550 -1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 16800 3282 168 168 16800 3282 16650 3207 -1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 18000 2400 150 150 18000 2400 18150 2400 -1 3 0 1 -1 -1 4 0 62 0.000 1 0.0000 16670 2400 168 168 16670 2400 16838 2568 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 1125 2400 150 150 1125 2400 1275 2550 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 4950 2400 150 150 4950 2400 5100 2550 +1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 5850 2400 150 150 5850 2400 6000 2400 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 4050 2400 150 150 4050 2400 4200 2550 +1 3 0 1 -1 7 0 0 0 0.000 1 0.0000 9375 1200 150 150 9375 1200 9525 1350 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 9375 3600 150 150 9375 3600 9525 3600 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 10050 2400 150 150 10050 2400 10200 2400 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 13725 2400 150 150 13725 2400 13875 2400 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 12375 2400 150 150 12375 2400 12525 2400 +1 3 2 1 -1 -1 0 0 -1 4.000 1 0.0000 9375 2400 75 75 9375 2400 9450 2475 +1 3 0 1 -1 7 0 0 0 0.000 1 0.0000 6600 1200 150 150 6600 1200 6750 1350 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 6600 3600 150 150 6600 3600 6750 3600 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 7275 2400 150 150 7275 2400 7425 2400 +1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 8175 2400 168 168 8175 2400 8325 2475 +1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 10950 2400 150 150 10950 2400 11100 2400 +1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 15300 1200 150 150 15300 1200 15450 1200 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 14400 1200 150 150 14400 1200 14550 1200 +1 3 0 1 -1 -1 0 0 20 0.000 1 0.0000 16290 2355 150 150 16290 2355 16440 2355 +1 3 0 1 7 -1 2 0 20 0.000 1 0.0000 17715 2730 150 150 17715 2730 17865 2880 +1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 17940 3180 168 168 17940 3180 18090 3255 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 17715 1980 150 150 17715 1980 17865 2130 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 19350 1593 168 168 19350 1593 19182 1425 +1 3 0 1 -1 -1 2 0 62 0.000 1 0.0000 19650 2400 150 150 19650 2400 19800 2550 +1 3 0 1 -1 -1 0 0 62 0.000 1 0.0000 19350 3282 168 168 19350 3282 19200 3207 +1 3 0 0 -1 -1 0 0 10 0.000 1 0.0000 20550 2400 150 150 20550 2400 20700 2400 +1 3 0 1 -1 -1 4 0 62 0.000 1 0.0000 19220 2400 168 168 19220 2400 19388 2568 +2 1 2 1 -1 7 1 0 -1 3.000 0 0 -1 0 0 2 + 750 2400 3300 2400 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2 + 4050 2400 4950 2400 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2 + 4950 2400 5850 2400 2 1 0 1 -1 -1 0 0 -1 0.000 0 0 -1 0 0 3 - 6825 1200 7500 2400 6825 3600 + 9375 1200 10050 2400 9375 3600 2 1 2 1 -1 -1 0 0 -1 3.000 0 0 -1 0 0 2 - 6825 900 6825 3900 + 9375 900 9375 3900 2 1 2 1 -1 -1 0 0 -1 3.000 0 0 -1 0 0 2 - 7500 2400 7500 1950 + 10050 2400 10050 1950 2 1 0 1 -1 7 1 0 -1 0.000 0 0 -1 0 0 2 - 7500 2400 8400 2400 + 10050 2400 10950 2400 2 1 2 1 -1 -1 1 0 -1 3.000 0 0 -1 0 0 2 - 8400 2400 8400 1950 + 10950 2400 10950 1950 2 1 2 1 -1 -1 1 0 -1 3.000 0 0 -1 0 0 2 - 6825 2400 8700 2400 + 9375 2400 11250 2400 2 1 2 1 -1 -1 0 0 -1 3.000 0 0 -1 0 0 2 - 7500 2100 8400 2100 + 10050 2100 10950 2100 2 1 2 1 -1 -1 0 0 -1 3.000 0 0 -1 0 0 2 - 11850 1200 11850 1650 + 14400 1200 14400 1650 2 1 0 1 -1 -1 1 0 -1 0.000 0 0 -1 0 0 4 - 9825 2400 11175 2400 11850 1200 12750 1200 + 12375 2400 13725 2400 14400 1200 15300 1200 2 1 2 1 -1 -1 1 0 -1 3.000 0 0 -1 0 0 2 - 12750 1200 12750 1650 + 15300 1200 15300 1650 2 1 2 1 -1 -1 0 0 -1 3.000 0 0 -1 0 0 2 - 11850 1500 12750 1500 + 14400 1500 15300 1500 2 1 2 1 -1 7 0 0 -1 3.000 0 0 -1 1 1 3 2 1 1.00 60.00 120.00 2 1 1.00 60.00 120.00 - 5175 1575 4725 2400 5175 3225 + 7725 1575 7275 2400 7725 3225 2 1 2 1 -1 7 1 0 -1 3.000 0 0 -1 0 0 3 - 5175 1575 5625 2400 5175 3225 + 7725 1575 8175 2400 7725 3225 2 1 0 1 -1 -1 0 0 -1 0.000 0 0 -1 0 0 3 - 4050 1200 4725 2400 4050 3600 + 6600 1200 7275 2400 6600 3600 2 1 0 1 -1 7 1 0 -1 0.000 0 0 -1 0 0 2 - 4725 2400 5625 2400 -2 1 2 1 -1 7 1 0 -1 3.000 0 0 -1 0 0 2 - 750 2400 3300 2400 + 7275 2400 8175 2400 2 1 0 2 -1 -1 1 0 -1 0.000 0 0 -1 0 1 2 1 1 1.00 70.00 660.00 - 15090 2355 15390 3180 + 17640 2355 17940 3180 2 3 0 1 7 -1 1 0 20 0.000 0 0 -1 0 0 4 - 15315 3180 15090 2355 15465 3180 15315 3180 + 17865 3180 17640 2355 18015 3180 17865 3180 2 3 2 1 -1 7 2 0 -1 3.000 0 0 -1 0 0 4 - 13740 2355 15165 1980 15165 2730 13740 2355 + 16290 2355 17715 1980 17715 2730 16290 2355 2 3 0 1 7 -1 1 0 20 0.000 0 0 -1 0 0 4 - 16870 3267 17095 2442 16720 3267 16870 3267 + 19420 3267 19645 2442 19270 3267 19420 3267 2 3 0 1 7 -1 1 0 20 0.000 0 0 -1 0 0 4 - 16875 1575 17100 2400 16725 1575 16875 1575 + 19425 1575 19650 2400 19275 1575 19425 1575 2 1 0 1 -1 7 1 0 -1 0.000 0 0 -1 0 0 2 - 17100 2400 18000 2400 + 19650 2400 20550 2400 2 1 2 1 -1 -1 0 0 -1 3.000 0 0 -1 0 0 2 - 17100 1950 17100 2400 + 19650 1950 19650 2400 2 1 2 1 -1 -1 1 0 -1 3.000 0 0 -1 0 0 2 - 18000 1950 18000 2400 + 20550 1950 20550 2400 2 1 2 1 -1 -1 0 0 -1 3.000 0 0 -1 0 0 2 - 17100 2100 18000 2100 + 19650 2100 20550 2100 2 3 0 0 0 7 3 0 50 0.000 0 0 -1 0 0 5 - 17025 2287 16650 2362 16650 2437 17025 2512 17025 2287 -4 0 -1 0 0 18 30 0.0000 4 315 255 1875 3600 2\001 -4 0 -1 0 0 17 24 0.0000 4 210 195 1522 2175 a\001 -4 0 -1 0 0 17 24 0.0000 4 270 615 2197 2175 1-a\001 -4 0 -1 0 0 17 24 0.0000 4 210 195 4747 1950 a\001 -4 0 -1 0 0 17 24 0.0000 4 285 210 4747 3150 b\001 -4 0 -1 0 0 17 24 0.0000 4 210 195 6547 1950 a\001 -4 0 -1 0 0 17 24 0.0000 4 270 615 6172 3075 1-a\001 -4 0 -1 0 0 18 30 0.0000 4 330 930 10725 3600 3fad\001 -4 0 -1 0 0 18 30 0.0000 4 315 255 4575 3600 3\001 -4 0 0 1 0 18 30 0.0000 4 330 675 7575 3600 3fd\001 -4 0 -1 0 0 32 24 0.0000 4 255 180 10950 1050 q\001 -4 0 -1 0 0 17 24 0.0000 4 285 210 12172 1875 d\001 -4 0 -1 1 0 16 24 0.0000 4 360 450 12045 1868 | |\001 -4 0 -1 0 0 17 24 0.0000 4 285 210 7818 1950 b\001 -4 0 -1 1 0 16 24 0.0000 4 360 450 7713 1950 | |\001 -4 0 -1 0 0 18 30 0.0000 4 315 945 13965 3555 3out\001 -4 0 -1 0 0 18 30 0.0000 4 330 675 17175 3600 4fd\001 -4 0 -1 0 0 17 24 0.0000 4 210 180 17422 2025 c\001 -4 0 -1 1 0 16 24 0.0000 4 360 450 17309 2025 | |\001 + 19575 2287 19200 2362 19200 2437 19575 2512 19575 2287 +4 0 -1 0 0 18 30 0.0000 4 345 270 1875 3600 2\001 +4 0 -1 0 0 17 24 0.0000 4 210 225 1522 2175 a\001 +4 0 -1 0 0 17 24 0.0000 4 285 585 2197 2175 1-a\001 +4 0 -1 0 0 17 24 0.0000 4 210 225 7297 1950 a\001 +4 0 -1 0 0 17 24 0.0000 4 300 225 7297 3150 b\001 +4 0 -1 0 0 17 24 0.0000 4 210 225 9097 1950 a\001 +4 0 -1 0 0 17 24 0.0000 4 285 585 8722 3075 1-a\001 +4 0 -1 0 0 18 30 0.0000 4 360 1005 13275 3600 3fad\001 +4 0 -1 0 0 18 30 0.0000 4 345 270 7125 3600 3\001 +4 0 0 1 0 18 30 0.0000 4 360 735 10125 3600 3fd\001 +4 0 -1 0 0 32 24 0.0000 4 285 210 13500 1050 q\001 +4 0 -1 0 0 17 24 0.0000 4 300 225 14722 1875 d\001 +4 0 -1 1 0 16 24 0.0000 4 390 570 14595 1868 | |\001 +4 0 -1 0 0 17 24 0.0000 4 300 225 10368 1950 b\001 +4 0 -1 1 0 16 24 0.0000 4 390 570 10263 1950 | |\001 +4 0 -1 0 0 18 30 0.0000 4 360 1035 16515 3555 3out\001 +4 0 -1 0 0 18 30 0.0000 4 360 735 19725 3600 4fd\001 +4 0 -1 0 0 17 24 0.0000 4 210 210 19972 2025 c\001 +4 0 -1 1 0 16 24 0.0000 4 390 570 19859 2025 | |\001 +4 0 -1 1 0 16 24 0.0000 4 390 570 5175 1950 | |\001 +4 0 -1 0 0 17 24 0.0000 4 210 225 5325 1950 a\001 +4 0 -1 0 0 18 30 0.0000 4 360 735 4575 3600 2fd\001 diff --git a/docs/reference-manual/functions/plots/dummies.pdf b/docs/reference-manual/functions/plots/dummies.pdf index 7835f7695d..1c8791a5d2 100644 Binary files a/docs/reference-manual/functions/plots/dummies.pdf and b/docs/reference-manual/functions/plots/dummies.pdf differ diff --git a/docs/reference-manual/functions/plots/dummies.svg b/docs/reference-manual/functions/plots/dummies.svg index 45fa90e44b..b009d62dd2 100644 --- a/docs/reference-manual/functions/plots/dummies.svg +++ b/docs/reference-manual/functions/plots/dummies.svg @@ -1,5344 +1,410 @@ - - - -image/svg+xml| | -3fd -| || | -1-aaba1-aa -23fad -3out4fdn -c -b -3 -θ -d - \ No newline at end of file + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +3fd + +| | + +| | + +| | + +| | + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +2 + +a + +1-a + +a + +b + +a + +1-a + +3fad + +3 + +θ + +d + +b + +3out + +4fd + +c + +a + +2fd + + diff --git a/docs/reference-manual/topologies/topology-file-formats.rst b/docs/reference-manual/topologies/topology-file-formats.rst index 084a585151..49401c2405 100644 --- a/docs/reference-manual/topologies/topology-file-formats.rst +++ b/docs/reference-manual/topologies/topology-file-formats.rst @@ -253,6 +253,8 @@ interactions can be converted to constraints by :ref:`grompp `. +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+ | 2-body virtual site | ``virtual_sites2`` | 3 | 1 | |AO| () | | +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+ + | 2-body virtual site (fd) | ``virtual_sites2`` | 3 | 2 | |DO| (nm) | | + +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+ | 3-body virtual site | ``virtual_sites3`` | 4 | 1 | |AO|, |BO| () | | +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+ | 3-body virtual site (fd) | ``virtual_sites3`` | 4 | 2 | |AO| (); |DO| (nm) | | diff --git a/docs/release-notes/2020/major/features.rst b/docs/release-notes/2020/major/features.rst index ed3ae87464..e27510b45c 100644 --- a/docs/release-notes/2020/major/features.rst +++ b/docs/release-notes/2020/major/features.rst @@ -17,3 +17,10 @@ increasing the similarity of a simulated density to the reference density. Multiple protocols are available for how to calculate simulated densities as well as how the similarity between a reference and a simulated density is evaluated. + +Virtual site on the line through two atoms at fixed distance +"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +This is use useful for e.g. halogens in the CHARMM force field. + +:issue:`2451` diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index c58d374902..300f56ba67 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -125,6 +125,7 @@ enum tpxv tpxv_MimicQMMM, /**< Introduced support for MiMiC QM/MM interface */ tpxv_PullAverage, /**< Added possibility to output average pull force and position */ tpxv_GenericInternalParameters, /**< Added internal parameters for mdrun modules*/ + tpxv_VSite2FD, /**< Added 2FD type virtual site */ tpxv_Count /**< the total number of tpxv versions */ }; @@ -207,6 +208,7 @@ static const t_ftupd ftupd[] = { { 79, F_DVDL_RESTRAINT }, { 79, F_DVDL_TEMPERATURE }, { tpxv_GenericInternalParameters, F_DENSITYFITTING }, + { tpxv_VSite2FD, F_VSITE2FD }, }; #define NFTUPD asize(ftupd) @@ -1968,6 +1970,7 @@ static void do_iparams(gmx::ISerializer *serializer, serializer->doReal(&iparams->settle.dhh); break; case F_VSITE2: + case F_VSITE2FD: serializer->doReal(&iparams->vsite.a); break; case F_VSITE3: diff --git a/src/gromacs/gmxlib/nrnb.cpp b/src/gromacs/gmxlib/nrnb.cpp index 08de8f3af1..1eea461821 100644 --- a/src/gromacs/gmxlib/nrnb.cpp +++ b/src/gromacs/gmxlib/nrnb.cpp @@ -174,6 +174,7 @@ static const t_nrnb_data nbdata[eNRNB] = { { "Constraint-Vir", 24 }, { "Settle", 323 }, { "Virtual Site 2", 23 }, + { "Virtual Site 2fd", 63 }, { "Virtual Site 3", 37 }, { "Virtual Site 3fd", 95 }, { "Virtual Site 3fad", 176 }, diff --git a/src/gromacs/gmxlib/nrnb.h b/src/gromacs/gmxlib/nrnb.h index 66decbebca..bab0ac594e 100644 --- a/src/gromacs/gmxlib/nrnb.h +++ b/src/gromacs/gmxlib/nrnb.h @@ -108,7 +108,8 @@ enum eNR_PCOUPL, eNR_EKIN, eNR_LINCS, eNR_LINCSMAT, eNR_SHAKE, eNR_CONSTR_V, eNR_SHAKE_RIJ, eNR_CONSTR_VIR, eNR_SETTLE, - eNR_VSITE2, eNR_VSITE3, eNR_VSITE3FD, + eNR_VSITE2, eNR_VSITE2FD, + eNR_VSITE3, eNR_VSITE3FD, eNR_VSITE3FAD, eNR_VSITE3OUT, eNR_VSITE4FD, eNR_VSITE4FDN, eNR_VSITEN, eNR_CMAP, eNR_UREY_BRADLEY, eNR_CROSS_BOND_BOND, diff --git a/src/gromacs/gmxpreprocess/convparm.cpp b/src/gromacs/gmxpreprocess/convparm.cpp index 16c77479cd..40cbc20ed1 100644 --- a/src/gromacs/gmxpreprocess/convparm.cpp +++ b/src/gromacs/gmxpreprocess/convparm.cpp @@ -406,6 +406,7 @@ assign_param(t_functype ftype, t_iparams *newparam, newparam->settle.dhh = old[1]; break; case F_VSITE2: + case F_VSITE2FD: case F_VSITE3: case F_VSITE3FD: case F_VSITE3OUT: diff --git a/src/gromacs/gmxpreprocess/topdirs.cpp b/src/gromacs/gmxpreprocess/topdirs.cpp index 89f4374b1f..8f7aea5bd7 100644 --- a/src/gromacs/gmxpreprocess/topdirs.cpp +++ b/src/gromacs/gmxpreprocess/topdirs.cpp @@ -207,7 +207,15 @@ int ifunc_index(Directive d, int type) return F_BHAM; } case Directive::d_vsites2: - return F_VSITE2; + switch (type) + { + case 1: + return F_VSITE2; + case 2: + return F_VSITE2FD; + default: + gmx_fatal(FARGS, "Invalid vsites2 type %d", type); + } case Directive::d_vsites3: switch (type) { diff --git a/src/gromacs/mdlib/vsite.cpp b/src/gromacs/mdlib/vsite.cpp index 063b63c9a3..cbc42a8f4b 100644 --- a/src/gromacs/mdlib/vsite.cpp +++ b/src/gromacs/mdlib/vsite.cpp @@ -189,6 +189,11 @@ static int pbc_rvec_sub(const t_pbc *pbc, const rvec xi, const rvec xj, rvec dx) } } +static inline real inverseNorm(const rvec x) +{ + return gmx::invsqrt(iprod(x, x)); +} + /* Vsite construction routines */ static void constr_vsite2(const rvec xi, const rvec xj, rvec x, real a, const t_pbc *pbc) @@ -215,6 +220,24 @@ static void constr_vsite2(const rvec xi, const rvec xj, rvec x, real a, const t_ /* TOTAL: 10 flops */ } +static void constr_vsite2FD(const rvec xi, const rvec xj, rvec x, real a, + const t_pbc *pbc) +{ + rvec xij; + pbc_rvec_sub(pbc, xj, xi, xij); + /* 3 flops */ + + const real b = a*inverseNorm(xij); + /* 6 + 10 flops */ + + x[XX] = xi[XX] + b*xij[XX]; + x[YY] = xi[YY] + b*xij[YY]; + x[ZZ] = xi[ZZ] + b*xij[ZZ]; + /* 6 Flops */ + + /* TOTAL: 25 flops */ +} + static void constr_vsite3(const rvec xi, const rvec xj, const rvec xk, rvec x, real a, real b, const t_pbc *pbc) { @@ -258,7 +281,7 @@ static void constr_vsite3FD(const rvec xi, const rvec xj, const rvec xk, rvec x, temp[ZZ] = xij[ZZ] + a*xjk[ZZ]; /* 6 flops */ - c = b*gmx::invsqrt(iprod(temp, temp)); + c = b*inverseNorm(temp); /* 6 + 10 flops */ x[XX] = xi[XX] + c*temp[XX]; @@ -278,13 +301,13 @@ static void constr_vsite3FAD(const rvec xi, const rvec xj, const rvec xk, rvec x pbc_rvec_sub(pbc, xk, xj, xjk); /* 6 flops */ - invdij = gmx::invsqrt(iprod(xij, xij)); + invdij = inverseNorm(xij); c1 = invdij * invdij * iprod(xij, xjk); xp[XX] = xjk[XX] - c1*xij[XX]; xp[YY] = xjk[YY] - c1*xij[YY]; xp[ZZ] = xjk[ZZ] - c1*xij[ZZ]; a1 = a*invdij; - b1 = b*gmx::invsqrt(iprod(xp, xp)); + b1 = b*inverseNorm(xp); /* 45 */ x[XX] = xi[XX] + a1*xij[XX] + b1*xp[XX]; @@ -330,7 +353,7 @@ static void constr_vsite4FD(const rvec xi, const rvec xj, const rvec xk, const r temp[ZZ] = xij[ZZ] + a*xjk[ZZ] + b*xjl[ZZ]; /* 12 flops */ - d = c*gmx::invsqrt(iprod(temp, temp)); + d = c*inverseNorm(temp); /* 6 + 10 flops */ x[XX] = xi[XX] + d*temp[XX]; @@ -369,7 +392,7 @@ static void constr_vsite4FDN(const rvec xi, const rvec xj, const rvec xk, const cprod(rja, rjb, rm); /* 9 flops */ - d = c*gmx::invsqrt(norm2(rm)); + d = c*inverseNorm(rm); /* 5+5+1 flops */ x[XX] = xi[XX] + d*rm[XX]; @@ -496,6 +519,10 @@ static void construct_vsites_thread(rvec x[], aj = ia[3]; constr_vsite2(x[ai], x[aj], x[avsite], a1, pbc_null2); break; + case F_VSITE2FD: + aj = ia[3]; + constr_vsite2FD(x[ai], x[aj], x[avsite], a1, pbc_null2); + break; case F_VSITE3: aj = ia[3]; ak = ia[4]; @@ -736,6 +763,105 @@ void constructVsitesGlobal(const gmx_mtop_t &mtop, } } +static void spread_vsite2FD(const t_iatom ia[], real a, + const rvec x[], + rvec f[], rvec fshift[], + gmx_bool VirCorr, matrix dxdf, + const t_pbc *pbc, const t_graph *g) +{ + const int av = ia[1]; + const int ai = ia[2]; + const int aj = ia[3]; + rvec fv; + copy_rvec(f[av], fv); + + rvec xij; + int sji = pbc_rvec_sub(pbc, x[aj], x[ai], xij); + /* 6 flops */ + + const real invDistance = inverseNorm(xij); + const real b = a*invDistance; + /* 4 + ?10? flops */ + + const real fproj = iprod(xij, fv)*invDistance*invDistance; + + rvec fj; + fj[XX] = b*(fv[XX] - fproj*xij[XX]); + fj[YY] = b*(fv[YY] - fproj*xij[YY]); + fj[ZZ] = b*(fv[ZZ] - fproj*xij[ZZ]); + /* 9 */ + + /* b is already calculated in constr_vsite2FD + storing b somewhere will save flops. */ + + f[ai][XX] += fv[XX] - fj[XX]; + f[ai][YY] += fv[YY] - fj[YY]; + f[ai][ZZ] += fv[ZZ] - fj[ZZ]; + f[aj][XX] += fj[XX]; + f[aj][YY] += fj[YY]; + f[aj][ZZ] += fj[ZZ]; + /* 9 Flops */ + + if (fshift) + { + int svi; + if (g) + { + ivec di; + ivec_sub(SHIFT_IVEC(g, ia[1]), SHIFT_IVEC(g, ai), di); + svi = IVEC2IS(di); + ivec_sub(SHIFT_IVEC(g, aj), SHIFT_IVEC(g, ai), di); + sji = IVEC2IS(di); + } + else if (pbc) + { + rvec xvi; + svi = pbc_rvec_sub(pbc, x[av], x[ai], xvi); + } + else + { + svi = CENTRAL; + } + + if (svi != CENTRAL || sji != CENTRAL) + { + rvec_dec(fshift[svi], fv); + fshift[CENTRAL][XX] += fv[XX] - fj[XX]; + fshift[CENTRAL][YY] += fv[YY] - fj[YY]; + fshift[CENTRAL][ZZ] += fv[ZZ] - fj[ZZ]; + fshift[ sji][XX] += fj[XX]; + fshift[ sji][YY] += fj[YY]; + fshift[ sji][ZZ] += fj[ZZ]; + } + } + + if (VirCorr) + { + /* When VirCorr=TRUE, the virial for the current forces is not + * calculated from the redistributed forces. This means that + * the effect of non-linear virtual site constructions on the virial + * needs to be added separately. This contribution can be calculated + * in many ways, but the simplest and cheapest way is to use + * the first constructing atom ai as a reference position in space: + * subtract (xv-xi)*fv and add (xj-xi)*fj. + */ + rvec xiv; + + pbc_rvec_sub(pbc, x[av], x[ai], xiv); + + for (int i = 0; i < DIM; i++) + { + for (int j = 0; j < DIM; j++) + { + /* As xix is a linear combination of j and k, use that here */ + dxdf[i][j] += -xiv[i]*fv[j] + xij[i]*fj[j]; + } + } + } + + /* TOTAL: 38 flops */ +} + static void spread_vsite3(const t_iatom ia[], real a, real b, const rvec x[], rvec f[], rvec fshift[], @@ -800,7 +926,7 @@ static void spread_vsite3FD(const t_iatom ia[], real a, real b, gmx_bool VirCorr, matrix dxdf, const t_pbc *pbc, const t_graph *g) { - real c, invl, fproj, a1; + real fproj, a1; rvec xvi, xij, xjk, xix, fv, temp; t_iatom av, ai, aj, ak; int svi, sji, skj; @@ -822,11 +948,11 @@ static void spread_vsite3FD(const t_iatom ia[], real a, real b, xix[ZZ] = xij[ZZ]+a*xjk[ZZ]; /* 6 flops */ - invl = gmx::invsqrt(iprod(xix, xix)); - c = b*invl; + const real invDistance = inverseNorm(xix); + const real c = b*invDistance; /* 4 + ?10? flops */ - fproj = iprod(xix, fv)*invl*invl; /* = (xix . f)/(xix . xix) */ + fproj = iprod(xix, fv)*invDistance*invDistance; /* = (xix . f)/(xix . xix) */ temp[XX] = c*(fv[XX]-fproj*xix[XX]); temp[YY] = c*(fv[YY]-fproj*xix[YY]); @@ -929,14 +1055,14 @@ static void spread_vsite3FAD(const t_iatom ia[], real a, real b, skj = pbc_rvec_sub(pbc, x[ak], x[aj], xjk); /* 6 flops */ - invdij = gmx::invsqrt(iprod(xij, xij)); + invdij = inverseNorm(xij); invdij2 = invdij * invdij; c1 = iprod(xij, xjk) * invdij2; xperp[XX] = xjk[XX] - c1*xij[XX]; xperp[YY] = xjk[YY] - c1*xij[YY]; xperp[ZZ] = xjk[ZZ] - c1*xij[ZZ]; /* xperp in plane ijk, perp. to ij */ - invdp = gmx::invsqrt(iprod(xperp, xperp)); + invdp = inverseNorm(xperp); a1 = a*invdij; b1 = b*invdp; /* 45 flops */ @@ -1122,7 +1248,7 @@ static void spread_vsite4FD(const t_iatom ia[], real a, real b, real c, gmx_bool VirCorr, matrix dxdf, const t_pbc *pbc, const t_graph *g) { - real d, invl, fproj, a1; + real fproj, a1; rvec xvi, xij, xjk, xjl, xix, fv, temp; int av, ai, aj, ak, al; ivec di; @@ -1146,13 +1272,13 @@ static void spread_vsite4FD(const t_iatom ia[], real a, real b, real c, } /* 12 flops */ - invl = gmx::invsqrt(iprod(xix, xix)); - d = c*invl; + const real invDistance = inverseNorm(xix); + const real d = c*invDistance; /* 4 + ?10? flops */ copy_rvec(f[av], fv); - fproj = iprod(xix, fv)*invl*invl; /* = (xix . f)/(xix . xix) */ + fproj = iprod(xix, fv)*invDistance*invDistance; /* = (xix . f)/(xix . xix) */ for (m = 0; m < DIM; m++) { @@ -1272,7 +1398,7 @@ static void spread_vsite4FDN(const t_iatom ia[], real a, real b, real c, cprod(rja, rjb, rm); /* 9 flops */ - invrm = gmx::invsqrt(norm2(rm)); + invrm = inverseNorm(rm); denom = invrm*invrm; /* 5+5+2 flops */ @@ -1482,6 +1608,9 @@ static void spread_vsite_f_thread(const rvec x[], case F_VSITE2: spread_vsite2(ia, a1, x, f, fshift, pbc_null2, g); break; + case F_VSITE2FD: + spread_vsite2FD(ia, a1, x, f, fshift, VirCorr, dxdf, pbc_null2, g); + break; case F_VSITE3: b1 = ip[tp].vsite.b; spread_vsite3(ia, a1, b1, x, f, fshift, pbc_null2, g); @@ -1740,6 +1869,7 @@ void spread_vsite_f(const gmx_vsite_t *vsite, } inc_nrnb(nrnb, eNR_VSITE2, vsite_count(idef->il, F_VSITE2)); + inc_nrnb(nrnb, eNR_VSITE2FD, vsite_count(idef->il, F_VSITE2FD)); inc_nrnb(nrnb, eNR_VSITE3, vsite_count(idef->il, F_VSITE3)); inc_nrnb(nrnb, eNR_VSITE3FD, vsite_count(idef->il, F_VSITE3FD)); inc_nrnb(nrnb, eNR_VSITE3FAD, vsite_count(idef->il, F_VSITE3FAD)); diff --git a/src/gromacs/topology/ifunc.cpp b/src/gromacs/topology/ifunc.cpp index 56346e70e9..ee60db622a 100644 --- a/src/gromacs/topology/ifunc.cpp +++ b/src/gromacs/topology/ifunc.cpp @@ -147,6 +147,7 @@ const t_interaction_function interaction_function[F_NRE] = def_shk ("CONSTRNC", "Constr. No Conn.", 2, 1, 1 ), def_shkcb ("SETTLE", "Settle", 3, 2, 0 ), def_vsite ("VSITE2", "Virtual site 2", 3, 1 ), + def_vsite ("VSITE2FD", "Virtual site 2fd", 3, 1 ), def_vsite ("VSITE3", "Virtual site 3", 4, 2 ), def_vsite ("VSITE3FD", "Virtual site 3fd", 4, 2 ), def_vsite ("VSITE3FAD", "Virtual site 3fad", 4, 2 ), diff --git a/src/gromacs/topology/ifunc.h b/src/gromacs/topology/ifunc.h index c3c1fffa43..4d4e5d3479 100644 --- a/src/gromacs/topology/ifunc.h +++ b/src/gromacs/topology/ifunc.h @@ -186,6 +186,7 @@ enum F_CONSTRNC, F_SETTLE, F_VSITE2, + F_VSITE2FD, F_VSITE3, F_VSITE3FD, F_VSITE3FAD,