2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
36 #include <types/simple.h>
38 #include "gmx_lapack.h"
41 F77_FUNC(slasq4,SLASQ4)(int *i0,
61 float gam, gap1, gap2;
72 nn = (*n0 << 2) + *pp;
75 if ( fabs(*dmin__ - *dn)<GMX_FLOAT_EPS*fabs(*dmin__ + *dn) ||
76 fabs(*dmin__ - *dn1)<GMX_FLOAT_EPS*fabs(*dmin__ + *dn1)) {
78 b1 = sqrt(z__[nn - 3]) * sqrt(z__[nn - 5]);
79 b2 = sqrt(z__[nn - 7]) * sqrt(z__[nn - 9]);
80 a2 = z__[nn - 7] + z__[nn - 5];
82 if ( fabs(*dmin__ - *dn)<GMX_FLOAT_EPS*fabs(*dmin__ + *dn) &&
83 fabs(*dmin1 - *dn1)<GMX_FLOAT_EPS*fabs(*dmin1 + *dn1)) {
85 gap2 = *dmin2 - a2 - *dmin2 * .25;
86 if (gap2 > 0. && gap2 > b2) {
87 gap1 = a2 - *dn - b2 / gap2 * b2;
89 gap1 = a2 - *dn - (b1 + b2);
91 if (gap1 > 0. && gap1 > b1) {
92 d__1 = *dn - b1 / gap1 * b1, d__2 = *dmin__ * .5;
93 s = (d__1>d__2) ? d__1 : d__2;
101 d__1 = s, d__2 = a2 - (b1 + b2);
102 s = (d__1<d__2) ? d__1 : d__2;
104 d__1 = s, d__2 = *dmin__ * .333;
105 s = (d__1>d__2) ? d__1 : d__2;
113 if (fabs(*dmin__ - *dn)<GMX_FLOAT_EPS*fabs(*dmin__ + *dn)) {
116 if (z__[nn - 5] > z__[nn - 7]) {
119 b2 = z__[nn - 5] / z__[nn - 7];
122 np = nn - (*pp << 1);
125 if (z__[np - 4] > z__[np - 2]) {
128 a2 = z__[np - 4] / z__[np - 2];
129 if (z__[nn - 9] > z__[nn - 11]) {
132 b2 = z__[nn - 9] / z__[nn - 11];
138 i__1 = (*i0 << 2) - 1 + *pp;
139 for (i4 = np; i4 >= i__1; i4 += -4) {
140 if (fabs(b2)<GMX_FLOAT_MIN) {
144 if (z__[i4] > z__[i4 - 2]) {
147 b2 *= z__[i4] / z__[i4 - 2];
149 if (((b2>b1) ? b2 : b1) * 100. < a2 || .563 < a2) {
158 s = gam * (1. - sqrt(a2)) / (a2 + 1.);
161 } else if (fabs(*dmin__ - *dn2)<GMX_FLOAT_EPS*fabs(*dmin__ + *dn2)) {
166 np = nn - (*pp << 1);
170 if (z__[np - 8] > b2 || z__[np - 4] > b1) {
173 a2 = z__[np - 8] / b2 * (z__[np - 4] / b1 + 1.);
177 b2 = z__[nn - 13] / z__[nn - 15];
179 i__1 = (*i0 << 2) - 1 + *pp;
180 for (i4 = nn - 17; i4 >= i__1; i4 += -4) {
181 if (fabs(b2)<GMX_FLOAT_MIN) {
185 if (z__[i4] > z__[i4 - 2]) {
188 b2 *= z__[i4] / z__[i4 - 2];
190 if (((b2>b1) ? b2 : b1) * 100. < a2 || .563 < a2) {
199 s = gam * (1. - sqrt(a2)) / (a2 + 1.);
204 g += (1. - g) * .333;
205 } else if (*ttype == -18) {
206 g = .083250000000000005;
214 } else if (*n0in == *n0 + 1) {
216 if ( fabs(*dmin1 - *dn1)<GMX_FLOAT_EPS*fabs(*dmin1 + *dn1) &&
217 fabs(*dmin2 - *dn2)<GMX_FLOAT_EPS*fabs(*dmin2 + *dn2)) {
221 if (z__[nn - 5] > z__[nn - 7]) {
224 b1 = z__[nn - 5] / z__[nn - 7];
226 if (fabs(b2)<GMX_FLOAT_MIN) {
229 i__1 = (*i0 << 2) - 1 + *pp;
230 for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
232 if (z__[i4] > z__[i4 - 2]) {
235 b1 *= z__[i4] / z__[i4 - 2];
237 if (((a2>b1) ? a2 : b1) * 100. < b2) {
242 b2 = sqrt(b2 * 1.05);
244 a2 = *dmin1 / (d__1 * d__1 + 1.);
245 gap2 = *dmin2 * .5 - a2;
246 if (gap2 > 0. && gap2 > b2 * a2) {
247 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
248 s = (d__1>d__2) ? d__1 : d__2;
250 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
251 s = (d__1>d__2) ? d__1 : d__2;
257 if (fabs(*dmin1 - *dn1)<GMX_FLOAT_EPS*fabs(*dmin1 + *dn1)) {
263 } else if (*n0in == *n0 + 2) {
265 if (fabs(*dmin2 - *dn2)<GMX_FLOAT_EPS*fabs(*dmin2 + *dn2) &&
266 z__[nn - 5] * 2. < z__[nn - 7]) {
269 if (z__[nn - 5] > z__[nn - 7]) {
272 b1 = z__[nn - 5] / z__[nn - 7];
274 if (fabs(b2)<GMX_FLOAT_MIN) {
277 i__1 = (*i0 << 2) - 1 + *pp;
278 for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
279 if (z__[i4] > z__[i4 - 2]) {
282 b1 *= z__[i4] / z__[i4 - 2];
284 if (b1 * 100. < b2) {
289 b2 = sqrt(b2 * 1.05);
291 a2 = *dmin2 / (d__1 * d__1 + 1.);
292 gap2 = z__[nn - 7] + z__[nn - 9] - sqrt(z__[nn - 11]) * sqrt(z__[
294 if (gap2 > 0. && gap2 > b2 * a2) {
295 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
296 s = (d__1>d__2) ? d__1 : d__2;
298 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
299 s = (d__1>d__2) ? d__1 : d__2;
305 } else if (*n0in > *n0 + 2) {