nb_free_energy.c

最著名最快的分子模拟软件

/*
 * $Id: nb_free_energy.c,v 1.7.2.3 2008/02/29 07:02:44 spoel Exp $
 * 
 *                This source code is part of
 * 
 *                 G   R   O   M   A   C   S
 * 
 *          GROningen MAchine for Chemical Simulations
 * 
 *                        VERSION 3.3.3
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2008, The GROMACS development team,
 * check out http://www.gromacs.org for more information.

 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * If you want to redistribute modifications, please consider that
 * scientific software is very special. Version control is crucial -
 * bugs must be traceable. We will be happy to consider code for
 * inclusion in the official distribution, but derived work must not
 * be called official GROMACS. Details are found in the README & COPYING
 * files - if they are missing, get the official version at www.gromacs.org.
 * 
 * To help us fund GROMACS development, we humbly ask that you cite
 * the papers on the package - you can find them in the top README file.
 * 
 * For more info, check our website at http://www.gromacs.org
 * 
 * And Hey:
 * Groningen Machine for Chemical Simulation
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <math.h>

#include "vec.h"
#include "typedefs.h"

void
gmx_nb_free_energy_kernel(int                  icoul,
                          int                  ivdw,
                          int                  nri,
                          int *                iinr,
                          int *                jindex,
                          int *                jjnr,
                          int *                shift,
                          real *               shiftvec,
                          real *               fshift,
                          int *                gid,
                          real *               x,
                          real *               f,
                          real *               chargeA,
                          real *               chargeB,
                          real                 facel,
                          real                 krf,
                          real                 crf,
                          real                 ewc,
                          real *               Vc,
                          int *                typeA,
                          int *                typeB,
                          int                  ntype,
                          real *               nbfp,
                          real *               Vvdw,
                          real                 tabscale,
                          real *               VFtab,
                          real                 lambda,
                          real *               dvdlambda,
                          real                 alpha,
			  int                  lam_power,
                          real                 def_sigma6,
                          int *                outeriter,
                          int *                inneriter)
{
    int           n,ii,is3,ii3,k,nj0,nj1,jnr,j3,ggid;
    real          shX,shY,shZ;
    real          Fscal,FscalA,FscalB,tx,ty,tz;
    real          VcoulA,VcoulB,VvdwA,VvdwB;
    real          rinv6,r,rt;
    real          iqA,iqB;
    real          qqA,qqB,vcoul,vctot,krsq;
    int           ntiA,ntiB;
    int           tjA,tjB;
    real          rinvsix;
    real          Vvdw6,Vvdwtot;
    real          Vvdw12;
    real          ix,iy,iz,fix,fiy,fiz;
    real          dx,dy,dz,rsq,r4,r6,rinv;
    real          c6A,c12A,c6B,c12B;
    real          dvdl,L1,alfA,alfB,dalfA,dalfB;
    real          sigma6a,sigma6b;
    real          rA,rinvA,rinv4A,rB,rinvB,rinv4B;
    int           do_coultab,do_vdwtab,do_tab,tab_elemsize;
    int           n0,n1,nnn;
    real          Y,F,G,H,Fp,Geps,Heps2,eps,eps2,VV,FF;
    double        isp=0.564189583547756;


    /* fix compiler warnings */
    nj1 = 0;
    n1  = 0;
    eps = 0;
    eps2 = 0;
   
    dvdl = 0;
    L1   = 1.0 - lambda;

    alfA  = alpha*(lam_power==2 ? lambda*lambda : lambda);
    alfB  = alpha*(lam_power==2 ? L1*L1 : L1);
    dalfA = alpha*lam_power/6.0*(lam_power==2 ? lambda : 1); 
    dalfB = alpha*lam_power/6.0*(lam_power==2 ? L1 : 1); 

    /* Ewald table is special (icoul==5) */
    
    do_coultab = (icoul==3);
    do_vdwtab  = (ivdw==3);
    
    do_tab = do_coultab || do_vdwtab;
    
    /* we always use the combined table here */
    tab_elemsize = 12;
    
    for(n=0; (n<nri); n++)
    {
        is3              = 3*shift[n];     
        shX              = shiftvec[is3];  
        shY              = shiftvec[is3+1];
        shZ              = shiftvec[is3+2];
        nj0              = jindex[n];      
        nj1              = jindex[n+1];    
        ii               = iinr[n];        
        ii3              = 3*ii;           
        ix               = shX + x[ii3+0];
        iy               = shY + x[ii3+1];
        iz               = shZ + x[ii3+2];
        iqA              = facel*chargeA[ii];
        iqB              = facel*chargeB[ii];
        ntiA             = 2*ntype*typeA[ii];
        ntiB             = 2*ntype*typeB[ii];
        vctot            = 0;              
        Vvdwtot          = 0;              
        fix              = 0;              
        fiy              = 0;              
        fiz              = 0;              
        
        for(k=nj0; (k<nj1); k++)
        {
            jnr              = jjnr[k];        
            j3               = 3*jnr;          
            dx               = ix - x[j3];      
            dy               = iy - x[j3+1];      
            dz               = iz - x[j3+2];      
            rsq              = dx*dx+dy*dy+dz*dz;
            rinv             = invsqrt(rsq);
            r                = rsq*rinv;
            tjA              = ntiA+2*typeA[jnr];
            tjB              = ntiB+2*typeB[jnr];
            c6A              = nbfp[tjA];
            c6B              = nbfp[tjB];
            c12A             = nbfp[tjA+1];
            c12B             = nbfp[tjB+1];
            qqA              = iqA*chargeA[jnr]; 
            qqB              = iqB*chargeB[jnr]; 
            
            if((c6A > 0) && (c12A > 0)) 
            {
                sigma6a           = c12A/c6A;
            }
            else 
            {
                sigma6a           = def_sigma6;
            }
            if((c6B > 0) && (c12B > 0))
            {
                sigma6b           = c12B/c6B;
            }
            else
            {
                sigma6b           = def_sigma6;
            }
                        
            r4               = rsq*rsq;
            r6               = r4*rsq;
            
            FscalA           = 0;
            VcoulA           = 0;
            VvdwA            = 0;
            rinv4A           = 0;
            
            /* Only spend time on A state if it is non-zero */
            if( (qqA != 0) || (c6A != 0) || (c12A != 0) ) 
            {
                rA             = pow(alfA*sigma6a+r6,1.0/6.0);
                rinvA          = 1.0/rA;
                rinv4A         = rinvA*rinvA;
                rinv4A         = rinv4A*rinv4A;

                
                if(do_tab)
                {
                    rt         = rA*tabscale;
                    n0         = rt;
                    eps        = rt-n0;
                    eps2       = eps*eps;
                    n1         = tab_elemsize*n0;
                }
                
                if(icoul==1 || icoul==5)
                {
                    /* simple cutoff */