My Documents/Codes/MD++.svn/trunk/src/ljbond.cpp

Go to the documentation of this file.

/*
  ljbond.cpp
  by Wei Cai  caiwei@mit.edu
  Last Modified : Fri Mar  2 18:51:39 2007

  FUNCTION  :  MD simulation package of Lennard-Jones (LJ) potential
               LJ parameters depend on species
               Atoms can also form covalent bond with neighbors

  Test cases: scripts/work/ar/lj2-md.script, ljbond.tcl
*/

#include "ljbond.h"

void LJBONDFrame::initparser()
{
    int i, j;
    char s[100];
    MDPARALLELFrame::initparser();

    /* input */
    for(i=0;i<MAXSP;i++)
       for(j=i;j<MAXSP;j++)
       {
           sprintf(s,"C12_%d%d",i,j);
           bindvar(s,&(_ALJ[i][j]),DOUBLE);
           sprintf(s,"C6_%d%d",i,j);
           bindvar(s,&(_BLJ[i][j]),DOUBLE);
       }

    bindvar("BOND_R0",&BOND_R0,DOUBLE);
    bindvar("BOND_K", &BOND_K, DOUBLE);
    bindvar("Rcut",&LJ_RC,DOUBLE);
}

int LJBONDFrame::exec(char *name)
{
    if(MDPARALLELFrame::exec(name)==0) return 0;
    bindcommand(name,"initLJ",initLJ());
    bindcommand(name,"makelipids",makelipids());
    bindcommand(name,"linklipids",linklipids());
    return -1;
}

void LJBONDFrame::Alloc()
{
    int size;
    MDPARALLELFrame::Alloc();

    size = _NP*allocmultiple;
    Realloc(num_bonds, int,size);
    Realloc(bond_index,int,size*MAXNUMBOND);

    memset(num_bonds, 0,sizeof(int)*size);
    memset(bond_index,0,sizeof(int)*size*MAXNUMBOND);
}


void LJBONDFrame::initvars()
{
    int i, j;
    for(i=0;i<MAXSP;i++)
       for(j=i;j<MAXSP;j++)
       {
           _ALJ[i][j] = 4;
           _BLJ[i][j] = 4;
       }

    BOND_R0 = LJ_LENGTH;
    BOND_K  = 2;
    LJ_RC = 2.37343077641 * LJ_LENGTH; /* initial value, 4th neighbor */
    _RLIST=LJ_RC+1.1;
    _SKIN=_RLIST-LJ_RC;
    MDPARALLELFrame::initvars();
}

void LJBONDFrame::initLJ()
{
    int i, j;

    for(i=0;i<MAXSP;i++)
       for(j=i;j<MAXSP;j++)
       {
           ALJ[i][j]=_ALJ[i][j]*LJ_ENERGY*POW12(LJ_LENGTH);
           BLJ[i][j]=_BLJ[i][j]*LJ_ENERGY*POW6(LJ_LENGTH);
            Uc[i][j]=(ALJ[i][j]*POW6(1/LJ_RC)-BLJ[i][j])*POW6(1/LJ_RC);
           DUDRc[i][j]=-(12*ALJ[i][j]/POW6(LJ_RC)-6*BLJ[i][j])/POW6(LJ_RC)/LJ_RC;
       }
    /* symmetrize LJ parameter matrix */
    for(i=1;i<MAXSP;i++)
       for(j=0;j<i;j++)
       {
           _ALJ[i][j] = _ALJ[j][i];
           _BLJ[i][j] = _BLJ[j][i];
            ALJ[i][j] =  ALJ[j][i];
            BLJ[i][j] =  BLJ[j][i];
             Uc[i][j] =   Uc[j][i];
          DUDRc[i][j] =DUDRc[j][i];
       }

    _RLIST=LJ_RC+1.1;
    _SKIN=_RLIST-LJ_RC;

}

void LJBONDFrame::calcprop()
{
    MDPARALLELFrame::calcprop();
}

void LJBONDFrame::lennard_jones_bond()
{
    int i,j,ipt,jpt,k, bonded;
    double U,r,r2,ri6, R0,U0,K0,ri02,ri06;
    double ALJ_local,BLJ_local,Uc_local,DUDRc_local;   
    Vector3 sij, rij, fij;
    
    DUMP(HIG"Lennard Jones"NOR);
        
    refreshneighborlist();
    
    _EPOT=0;

    for(i=0;i<_NP;i++)
    {_F[i].clear(); _EPOT_IND[i]=0;}
    _VIRIAL.clear();

    for(ipt=0;ipt<_NP;ipt++)
    {
        for(j=0;j<nn[ipt];j++)
        {
            jpt=nindex[ipt][j];
            if(ipt>jpt) continue;

            /* test if atoms ipt and jpt have a covalent bond */
            bonded = 0;
            for(k=0;k<num_bonds[ipt];k++)
            {
                if(bond_index[ipt*MAXNUMBOND+k] == jpt)
                {
                    bonded = 1;
                    break;
                }
            }
            /* exclude covalently bonded pairs from LJ interaction */
            if(bonded) continue;

            sij=_SR[jpt]-_SR[ipt];
            sij.subint();
            rij=_H*sij;
            r2=rij.norm2();
            r=sqrt(r2);
            if(r<=LJ_RC)
            {
                ALJ_local =   ALJ[species[ipt]][species[jpt]];
                BLJ_local =   BLJ[species[ipt]][species[jpt]];
                 Uc_local =    Uc[species[ipt]][species[jpt]];
              DUDRc_local = DUDRc[species[ipt]][species[jpt]];

                R0 = 0.8;  /* short range cut-off */
                if(r<R0)
                {
                  ri02  = 1./(R0*R0);
                  ri06 = ri02*ri02*ri02;
                  U0 = (7*ALJ_local*ri06 - 4*BLJ_local)*ri06; 
                  K0 = (6*ALJ_local*ri06 - 3*BLJ_local)*ri06*ri02;
                  U = U0 - K0*r2 - r*DUDRc_local + (LJ_RC*DUDRc_local-Uc_local);
                  fij = rij*(2*K0 + DUDRc_local/r);
                }
                else {

                   ri6=1./(r2*r2*r2);
                
                   U=(ALJ_local*ri6-BLJ_local)*ri6-r*DUDRc_local+(LJ_RC*DUDRc_local-Uc_local);
                   fij=rij*((12*ALJ_local*ri6-6*BLJ_local)*ri6/r2+DUDRc_local/r);
                }

                _F[ipt]-=fij;
                _F[jpt]+=fij;
                _EPOT_IND[ipt]+=U*0.5;
                _EPOT_IND[jpt]+=U*0.5;
                _EPOT+=U;
                _VIRIAL.addnvv(1.,fij,rij);
            }
        }
    }

    /* covalently bonded interaction */
  if(BOND_K != 0)
  {
    for(ipt=0;ipt<_NP;ipt++)
    {
       //INFO_Printf("ljbond: ipt=%d num_bonds = %d\n",ipt,num_bonds[ipt]);
       for(j=0;j<num_bonds[ipt];j++)
       {
          jpt = bond_index[ipt*MAXNUMBOND + j];
          //INFO_Printf("ljbond: ipt=%d jpt=%d\n",ipt,jpt);

          if(ipt>jpt) continue;

          sij=_SR[jpt]-_SR[ipt];
          sij.subint();
          rij=_H*sij;
          r2=rij.norm2();
          r=sqrt(r2);
                
          U = BOND_K * SQR(r-BOND_R0);
          fij = rij * (-2*BOND_K*(1-BOND_R0/r)); /* bug fix from old ljbond.cpp */

          _F[ipt]-=fij;
          _F[jpt]+=fij;
          _EPOT_IND[ipt]+=U*0.5;
          _EPOT_IND[jpt]+=U*0.5;
          _EPOT+=U;
          _VIRIAL.addnvv(1.,fij,rij);
       }
       //INFO_Printf("ljbond: F[%d] = (%e %e %e)\n",ipt,_F[ipt].x,_F[ipt].y,_F[ipt].z);       
    }
  }
}

void LJBONDFrame::linklipids()
{
    int nsolvent, nlipid, chainlen, headid;
    int ipt, i, j;
    double bond_len;
    Vector3 dr, dr1, ds;
    Matrix33 hinv;

    nsolvent = (int) input[0];
    nlipid = (int) input[1];
    chainlen = (int) input[2];
    headid = (int) input[3];
    bond_len = input[4];

    ipt = 0;
    for(i=0;i<nlipid;i++)
    {
        for(j=0;j<chainlen;j++)
        {
           if(j==0)
           {
               /* only one bond */
               num_bonds[ipt] = 1;
               bond_index[ipt*MAXNUMBOND + 0] = ipt+1;
               
           } else if (j==(chainlen-1))
           {   /* only one bond */
               num_bonds[ipt]=1;
               bond_index[ipt*MAXNUMBOND + 0] = ipt-1;
           } else
           {   /* two bonds */
               num_bonds[ipt]=2;
               bond_index[ipt*MAXNUMBOND + 0] = ipt-1;
               bond_index[ipt*MAXNUMBOND + 1] = ipt+1;
           }
           ipt ++;
        }
    }
}

void LJBONDFrame::makelipids()
{
    int nsolvent, nlipid, chainlen, headid;
    int ipt, i, j;
    double bond_len;
    Vector3 dr, dr1, ds;
    Matrix33 hinv;

    nsolvent = (int) input[0];
    nlipid = (int) input[1];
    chainlen = (int) input[2];
    headid = (int) input[3];
    bond_len = input[4];

    _H.clear();
    _H[0][0] = input[5];
    _H[1][1] = input[6];
    _H[2][2] = input[7];

    _NP = nlipid*chainlen + nsolvent;
    Alloc();

    hinv = _H.inv();
    ipt = 0;
    for(i=0;i<nlipid;i++)
    {
        for(j=0;j<chainlen;j++)
        {
           if(j==0)
           {
               /* randomly position the first atom */
               _SR[ipt].set(drand48()-0.5,drand48()-0.5,drand48()-0.5);

               /* randomly choose a growth direction */
               dr.set(drand48()-0.5,drand48()-0.5,drand48()-0.5);
               dr/=dr.norm();  dr*=bond_len; ds = hinv*dr;

           } else {
               /* add a new atom in the randomly chosen direction */
               _SR[ipt] = _SR[ipt-1] + ds;
           }

           if(j==headid)
           {
               species[ipt] = 2; /* head particle (hydrophilic) */

               /* choose a slightly different growth direction */
               dr1.set(drand48()-0.5,drand48()-0.5,drand48()-0.5);
               dr1/=dr1.norm(); dr1*=(bond_len*0.5); dr*=-1; dr+=dr1;
               dr/=dr.norm();  dr*=bond_len; ds = hinv*dr;

           } else {
               species[ipt] = 1; /* tail particle (hydrophobic) */
           }

           ipt ++;
        }
    }
    linklipids();
     
    for(i=0;i<nsolvent;i++)
    {
        /* randomly position the solvent particle */
        _SR[ipt].set(drand48()-0.5,drand48()-0.5,drand48()-0.5);
        species[ipt] = 0;     /* solvent particle (water) */
        ipt ++;
    }             
}

void LJBONDFrame::potential()
{
    lennard_jones_bond();
}

#ifdef _TEST

/* Main Program Begins */
class LJBONDFrame sim;

/* The main program is defined here */
#include "main.cpp"

#endif//_TEST

---