nb_free_energy.c

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

                    VcoulA     = qqA*rinvA;
                    FscalA     = VcoulA*rinvA*rinvA;
                }
                else if(icoul==2)
                {
                    /* reaction-field */
                    krsq       = krf*rA*rA;      
                    VcoulA     = qqA*(rinvA+krsq-crf);
                    FscalA     = qqA*(rinvA-2.0*krsq)*rinvA*rinvA;
                }
                else if(icoul==3)
                {
                    /* non-Ewald tabulated coulomb */
                    nnn        = n1;
                    Y          = VFtab[nnn];
                    F          = VFtab[nnn+1];
                    Geps       = eps*VFtab[nnn+2];
                    Heps2      = eps2*VFtab[nnn+3];
                    Fp         = F+Geps+Heps2;
                    VV         = Y+eps*Fp;
                    FF         = Fp+Geps+2.0*Heps2;
                    VcoulA     = qqA*VV;
                    FscalA     = -qqA*tabscale*FF*rinvA;                    
                }
                
                if(ivdw==1)
                {
                    /* cutoff LJ */
                    rinv6            = rinvA*rinvA*rinv4A;
                    Vvdw6            = c6A*rinv6;     
                    Vvdw12           = c12A*rinv6*rinv6;
                    VvdwA            = Vvdw12-Vvdw6;
                    FscalA          += (12.0*Vvdw12-6.0*Vvdw6)*rinvA*rinvA;                    
                }
                else if(ivdw==3)
                {
                    /* Table LJ */
		    nnn = n1+4;
                    
                    /* dispersion */
                    Y          = VFtab[nnn];
                    F          = VFtab[nnn+1];
                    Geps       = eps*VFtab[nnn+2];
                    Heps2      = eps2*VFtab[nnn+3];
                    Fp         = F+Geps+Heps2;
                    VV         = Y+eps*Fp;
                    FF         = Fp+Geps+2.0*Heps2;
                    VvdwA     += c6A*VV;
                    FscalA    -= c6A*tabscale*FF*rinvA;                    
                    
                    /* repulsion */
                    Y          = VFtab[nnn+4];
                    F          = VFtab[nnn+5];
                    Geps       = eps*VFtab[nnn+6];
                    Heps2      = eps2*VFtab[nnn+7];
                    Fp         = F+Geps+Heps2;
                    VV         = Y+eps*Fp;
                    FF         = Fp+Geps+2.0*Heps2;
                    VvdwA     += c12A*VV;
                    FscalA    -= c12A*tabscale*FF*rinvA;
                }           
                /* Buckingham vdw free energy not supported */
            }
            
            FscalB           = 0;
            VcoulB           = 0;
            VvdwB            = 0;
            rinv4B           = 0;
            
            /* Only spend time on B state if it is non-zero */
            if( (qqB != 0) || (c6B != 0) || (c12B != 0) ) 
            {
                rB             = pow(alfB*sigma6b+r6,1.0/6.0);
                rinvB          = 1.0/rB;
                rinv4B         = rinvB*rinvB;
                rinv4B         = rinv4B*rinv4B;
                
                
                if(do_tab)
                {
                    rt         = rB*tabscale;
                    n0         = rt;
                    eps        = rt-n0;
                    eps2       = eps*eps;
                    n1         = tab_elemsize*n0;
                }
                
                if(icoul==1 || icoul==5)
                {
                    /* simple cutoff */
                    VcoulB     = qqB*rinvB;
                    FscalB     = VcoulB*rinvB*rinvB;
                }
                else if(icoul==2)
                {
                    /* reaction-field */
                    krsq       = krf*rB*rB;      
                    VcoulB     = qqB*(rinvB+krsq-crf);
                    FscalB     = qqB*(rinvB-2.0*krsq)*rinvB*rinvB;                    
                }
                else if(icoul==3)
                {
                    /* non-Ewald tabulated coulomb */
                    nnn        = n1;
                    Y          = VFtab[nnn];
                    F          = VFtab[nnn+1];
                    Geps       = eps*VFtab[nnn+2];
                    Heps2      = eps2*VFtab[nnn+3];
                    Fp         = F+Geps+Heps2;
                    VV         = Y+eps*Fp;
                    FF         = Fp+Geps+2.0*Heps2;
                    VcoulB     = qqB*VV;
                    FscalB     = -qqB*tabscale*FF*rinvB;                    
                }
                
                if(ivdw==1)
                {
                    /* cutoff LJ */
                    rinv6            = rinvB*rinvB*rinv4B;
                    Vvdw6            = c6B*rinv6;     
                    Vvdw12           = c12B*rinv6*rinv6;
                    VvdwB            = Vvdw12-Vvdw6;
                    FscalB          += (12.0*Vvdw12-6.0*Vvdw6)*rinvB*rinvB;                    
                }
                else if(ivdw==3)
                {
                    /* Table LJ */
                    nnn = n1+4;
                    
                    /* dispersion */
                    Y          = VFtab[nnn];
                    F          = VFtab[nnn+1];
                    Geps       = eps*VFtab[nnn+2];
                    Heps2      = eps2*VFtab[nnn+3];
                    Fp         = F+Geps+Heps2;
                    VV         = Y+eps*Fp;
                    FF         = Fp+Geps+2.0*Heps2;
                    VvdwB     += c6B*VV;
                    FscalB    -= c6B*tabscale*FF*rinvB;                    
                    
                    /* repulsion */
                    Y          = VFtab[nnn+4];
                    F          = VFtab[nnn+5];
                    Geps       = eps*VFtab[nnn+6];
                    Heps2      = eps2*VFtab[nnn+7];
                    Fp         = F+Geps+Heps2;
                    VV         = Y+eps*Fp;
                    FF         = Fp+Geps+2.0*Heps2;
                    VvdwB     += c12B*VV;
                    FscalB    -= c12B*tabscale*FF*rinvB;                    
                }           
                /* Buckingham vdw free energy not supported */
            }

            Fscal = 0;
            
            if(icoul==5)
            {
                /* Soft-core Ewald interactions are special:
                 * For the direct space interactions we effectively want the
                 * normal coulomb interaction (added above when icoul==5),
                 * but need to subtract the part added in reciprocal space.
                 */
                if (r != 0) 
                {
                    VV    = erf(ewc*r)*rinv;
                    FF    = rinv*rinv*(VV - 2.0*ewc*isp*exp(-ewc*ewc*rsq));
                }
                else 
                {
                    VV    = ewc*2.0/sqrt(M_PI);
                    FF    = 0;
                }
                vctot  -= (lambda*qqB + L1*qqA)*VV;
                Fscal  -= (lambda*qqB + L1*qqA)*FF;
                dvdl   -= (qqB - qqA)*VV;
            }
	    
            /* Assemble A and B states */
            vctot         += lambda*VcoulB + L1*VcoulA;
            Vvdwtot       += lambda*VvdwB  + L1*VvdwA;
                
            Fscal         += (L1*FscalA*rinv4A + lambda*FscalB*rinv4B)*r4;
            dvdl          += (VcoulB + VvdwB) - (VcoulA + VvdwA);
            dvdl          += lambda*dalfB*FscalB*sigma6b*rinv4B
	                       - L1*dalfA*FscalA*sigma6a*rinv4A;
                
            tx             = Fscal*dx;     
            ty             = Fscal*dy;     
            tz             = Fscal*dz;     
            fix            = fix + tx;      
            fiy            = fiy + ty;      
            fiz            = fiz + tz;      
            f[j3]          = f[j3]   - tx;
            f[j3+1]        = f[j3+1] - ty;
            f[j3+2]        = f[j3+2] - tz;
        }
        
        f[ii3]           = f[ii3]        + fix;
        f[ii3+1]         = f[ii3+1]      + fiy;
        f[ii3+2]         = f[ii3+2]      + fiz;
        fshift[is3]      = fshift[is3]   + fix;
        fshift[is3+1]    = fshift[is3+1] + fiy;
        fshift[is3+2]    = fshift[is3+2] + fiz;
        ggid             = gid[n];         
        Vc[ggid]         = Vc[ggid] + vctot;
        Vvdw[ggid]       = Vvdw[ggid] + Vvdwtot;
    }
    
    *dvdlambda      += dvdl;
    *outeriter       = nri;            
    *inneriter       = nj1;            
}