nonbonded.c

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

                
                gmx_nb_free_energy_kernel(nlist->icoul,
                                          nlist->ivdw,
                                          nlist->nri,
                                          nlist->iinr,
                                          nlist->jindex,
                                          nlist->jjnr,
                                          nlist->shift,
                                          fr->shift_vec[0],
                                          fshift,
                                          nlist->gid,
                                          x[0],
                                          f[0],
                                          mdatoms->chargeA,
                                          mdatoms->chargeB,
                                          fr->epsfac,
                                          fr->k_rf,
                                          fr->c_rf,
                                          fr->ewaldcoeff,
                                          egcoul,
                                          mdatoms->typeA,
                                          mdatoms->typeB,
                                          fr->ntype,
                                          fr->nbfp,
                                          egnb,
                                          nblists->tab.scale,
                                          tabledata,
                                          lambda,
                                          dvdlambda,
                                          fr->sc_alpha,
					  fr->sc_power,
                                          fr->sc_sigma6,
                                          &outeriter,
                                          &inneriter);
            }
            else
            {
	  /* Not free energy - call nonbonded kernel from function pointer */
	  kernelptr = nb_kernel_list[nrnb_ind];
	  
                if(kernelptr == NULL)
                {
                    /* Call generic nonbonded kernel */
                    /*
                    gmx_nb_generic_kernel(nlist,
                                          x,
                                          fr,
                                          f,
                                          fshift,
                                          mdatoms,
                                          lambda,
                                          dvdlambda,
                                          egcoul,
                                          egnb,
                                          &(fr->tab.scale),
                                          tabledata,
                                          &nthreads,
                                          &(nlist->count),
                                          nlist->mtx,
                                          &outeriter,
                                          &inneriter);
                     */
	  }
	  
	  /* Call the appropriate nonbonded kernel function */
	  (*kernelptr)( &(nlist->nri),
			nlist->iinr,
			nlist->jindex,
			nlist->jjnr,
			nlist->shift,
			fr->shift_vec[0],
			fshift,
			nlist->gid,
			x[0],
			f[0],
			mdatoms->chargeA,
			&(fr->epsfac),
			&(fr->k_rf),
			&(fr->c_rf),
			egcoul,
			mdatoms->typeA,
			&(fr->ntype),
			fr->nbfp,
			egnb,
			&(nblists->tab.scale),
			tabledata,
			NULL,
			NULL,
			NULL,
			NULL,
			&nthreads,
			&(nlist->count),
			nlist->mtx,
			&outeriter,
			&inneriter,
			NULL);
	}
	
	/* Update flop accounting */
	
	/* Outer loop in kernel */
	if (nlist->solvent_opt==enlistWATER)
            {
	  inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,wateratoms*outeriter);
            }
	else if (nlist->solvent_opt==enlistWATERWATER)
            {
	  inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,wateratoms*wateratoms*outeriter);
            }
	else
            {
	  inc_nrnb(nrnb,eNR_NBKERNEL_OUTER,outeriter);
            }
	/* inner loop in kernel */
	inc_nrnb(nrnb,nrnb_ind,inneriter);
      }
    }
  }
}


real 
do_nonbonded14(int nbonds,const t_iatom iatoms[],const t_iparams iparams[],
               const rvec x[],rvec f[],rvec fshift[],
               const t_pbc *pbc,const t_graph *g,
               real lambda,real *dvdlambda,
               const t_mdatoms *md,
               const t_forcerec *fr,int ngrp,real egnb[],real egcoul[])
{
    static    bool bWarn=FALSE;
    bool      bFullPBC;
    real      eps,r2,rtab2;
    rvec      dx,x14[2],f14[2];
    int       i,ai,aj,itype;
    int       typeA[2]={0,0},typeB[2]={0,1};
    real      chargeA[2],chargeB[2];
    int       gid,shift_vir,shift_f;
    int       j_index[] = { 0, 1 };
    int       i0=0,i1=1,i2=2;
    ivec      dt;
    int       outeriter,inneriter;
    int       nthreads = 1;
    int       count;
    real      krf,crf,tabscale;
    real *    nbfp;
    t_nblist  tmplist;
    int       icoul,ivdw;
    
#if GMX_THREADS
    pthread_mutex_t mtx;
#else
    void *    mtx = NULL;
#endif

    
#if GMX_THREADS
    pthread_mutex_initialize(&mtx);
#endif

    krf = fr->k_rf;
    crf = fr->c_rf;
    tabscale = fr->tab14.scale;
    
    /* Determine the values for icoul/ivdw. */
    if (fr->bEwald) {
      icoul = 1;
    } 
    else if(fr->bcoultab)
    {
        icoul = 3;
    }
    else if(fr->eeltype == eelRF_NEC)
    {
        icoul = 2;
    }
    else 
    {
        icoul = 1;
    }
    
    if(fr->bvdwtab)
    {
        ivdw = 3;
    }
    else if(fr->bBHAM)
    {
        ivdw = 2;
    }
    else 
    {
        ivdw = 1;
    }
    
    
    /* We don't do SSE or altivec here, due to large overhead for 4-fold 
     * unrolling on short lists 
     */
    
    bFullPBC = (fr->ePBC == epbcFULL);
    
    /* Reaction field stuff */  
    eps   = fr->epsfac*fr->fudgeQQ;
    
    rtab2 = sqr(fr->tab14.r);
    
    for(i=0; (i<nbonds); ) 
    {
        itype = iatoms[i++];
        ai    = iatoms[i++];
        aj    = iatoms[i++];

        if (!bFullPBC) 
        {
            /* This is a bonded interaction, atoms are in the same box */
            shift_f = CENTRAL;
            r2 = distance2(x[ai],x[aj]);
        }
        else 
        {
            /* Apply full periodic boundary conditions */
            shift_f = pbc_dx(pbc,x[ai],x[aj],dx);
            r2 = norm2(dx);
        }

        if (r2 >= rtab2) 
        {
            if (!bWarn) 
            {
                fprintf(stderr,"Warning: 1-4 interaction between %d and %d "
                        "at distance %.3f which is larger than the 1-4 table size %.3f nm\n", ai+1, aj+1, sqrt(r2), sqrt(rtab2));
                fprintf(stderr,"These are ignored for the rest of the simulation\n");
		fprintf(stderr,"This usually means your system is exploding,\n"
			"if not, you should increase table-extension in your mdp file\n");
                bWarn = TRUE;
            }
            if (debug) 
	      fprintf(debug,"%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored",
		      x[ai][XX],x[ai][YY],x[ai][ZZ],
		      x[aj][XX],x[aj][YY],x[aj][ZZ],
		      (int)ai+1,(int)aj+1,sqrt(r2));
        }
        else 
        {
            chargeA[0] = md->chargeA[ai];
            chargeA[1] = md->chargeA[aj];
            gid  = GID(md->cENER[ai],md->cENER[aj],ngrp);
            copy_rvec(x[ai],x14[0]);
            copy_rvec(x[aj],x14[1]);
            clear_rvec(f14[0]);
            clear_rvec(f14[1]);
#ifdef DEBUG
            fprintf(debug,"LJ14: grp-i=%2d, grp-j=%2d, ngrp=%2d, GID=%d\n",
                    md->cENER[ai],md->cENER[aj],ngrp,gid);
#endif
            
	    outeriter = inneriter = count = 0;
            nbfp = (real *)&(iparams[itype].lj14.c6A);

            if (fr->efep != efepNO &&
                (md->bPerturbed[ai] || md->bPerturbed[aj] ||
                 iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
                 iparams[itype].lj14.c12A != iparams[itype].lj14.c12B)) 
	      {
                
                chargeB[0] = md->chargeB[ai];
                chargeB[1] = md->chargeB[aj];
                /* We pass &(iparams[itype].lj14.c6A) as LJ parameter matrix
                 * to the innerloops.
                 * Here we use that the LJ-14 parameters are stored in iparams
                 * as c6A,c12A,c6B,c12B, which are referenced correctly
                 * in the innerloops if we assign type combinations 0-0 and 0-1
                 * to atom pair ai-aj in topologies A and B respectively.
                 */
                if(ivdw==2)
                {
                    gmx_fatal(FARGS,"Cannot do free energy Buckingham interactions.");
                }
                count = 0;
                gmx_nb_free_energy_kernel(icoul,
                                          ivdw,
                                          i1,
                                          &i0,
                                          j_index,
                                          &i1,
                                          &shift_f,
                                          fr->shift_vec[0],
                                          fshift[0],
                                          &gid,
                                          x14[0],
                                          f14[0],
                                          chargeA,
                                          chargeB,
                                          eps,
                                          krf,
                                          crf,
					  fr->ewaldcoeff,
                                          egcoul,
                                          typeA,
                                          typeB,
                                          i1,
                                          nbfp,
                                          egnb,
                                          tabscale,
                                          fr->tab14.tab,
                                          lambda,
                                          dvdlambda,
                                          fr->sc_alpha,
					  fr->sc_power,
                                          fr->sc_sigma6,
                                          &outeriter,
                                          &inneriter);
            }
            else 
            { 
                /* Not perturbed - call kernel 330 */
                F77_OR_C_FUNC_(nb_kernel330,NB_KERNEL330)
                    ( &i1,
                      &i0,
                      j_index,
                      &i1,
                      &shift_f,
                      fr->shift_vec[0],
                      fshift[0],
                      &gid,
                      x14[0],
                      f14[0],
                      chargeA,
                      &eps,
                      &krf,
                      &crf,
                      egcoul,
                      typeA,
                      &i1,
                      nbfp,
                      egnb,
                      &tabscale,
                      fr->tab14.tab,
                      NULL,
                      NULL,
                      NULL,
                      NULL,
                      &nthreads,
                      &count,
                      (void *)&mtx,
                      &outeriter,
                      &inneriter,
                      NULL);                
            }
            
            /* Add the forces */
            rvec_inc(f[ai],f14[0]);
            rvec_dec(f[aj],f14[0]);
            
            if (g) 
            {
                /* Correct the shift forces using the graph */
                ivec_sub(SHIFT_IVEC(g,ai),SHIFT_IVEC(g,aj),dt);    
                shift_vir = IVEC2IS(dt);
                rvec_inc(fshift[shift_vir],f14[0]);
                rvec_dec(fshift[CENTRAL],f14[0]);
            }

	    /* flops: eNR_KERNEL_OUTER + eNR_KERNEL330 + 12 */
        }
    }
    return 0.0;
}