// Last Modified : Wed Jun  4 21:13:52 1997
// $Id$


/*================================================
 Change Log::

Mon May  5 22:17:23 1997: 1.00
   First runnable version, almost(if not complete) bug free
Tue May  6 15:38:52 1997: 1.01
   Fixed a bug in NeibList, and FILEs
   Starting to make memory pool.
Wed May  7 00:08:19 1997: 1.02
   Delete all BC conditions for performance
Mon May 12 15:30:42 1997: 1.10
   Now works fine and fast for NeibList, start to add BinList
   Adding BinList, only support 2DIM pipe now.
Mon May 12 21:37:13 1997: 1.20
   BinList works, and NeibList works as well. Start to make shared mem
Mon May 12 23:00:33 1997: 1.30
   Shared memory works, start multiprocess [todo: wall(as Force Single)]
Wed May 28 18:41:26 1997: 1.4.0
   Get everything work and write good documents in `spec.h'
   Restrictions: only 2d,
                 only x binplane RPM,
                 only one block of wall region, must be rectangle, sc lattice
                 only divde processes horizontally along a pipe
   Minor bugs:first property dump is not right!
*/
#define _VERSION "1.4.0"
/*===============================================*/


#include <unistd.h>
#include <stdio.h>
#include <fcntl.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "Mng.h"


#define _S static
#define SQ(x) ((x)*(x))
//======================================================
// User Specifications
//======================================================
#include "spec.h"
//======================================================
// CPP Conditions
//======================================================

//We use macro `_3' to facilitate 2/3 dimesional compatible program
//Since a macro in C/C++ can only have definite number of parameters
//  we use `_' instead of `,' to make `_3' more useful.
#ifndef _2DIM
#define _V (_CUBEX*_CUBEY*_CUBEZ)
#define _DIMINT 3
#define _ ,
#define _3(x) x
#else
#define _V (_CUBEX*_CUBEY)
#define _DIMINT 2
#define _3(x) 
#endif

//_V: Volume
//_DR: Determined number density
#define _DR (_NP/_V)
      
//default value of _LCFG_X filenames
#ifndef _LCFG_READ
#define _LCFG_READ "config.dat"
#endif
#ifndef _LCFG_WRITE
#define _LCFG_WRITE "config.dat"
#endif

//default no rescale at all
#ifndef _NRESCALE
#define _NRESCALE 0
#endif


//-------- Neib part -----------

//define for optimized? flash distance
#ifdef _2DIM
#define _W 0.1
#else
#define _W 0.16
#endif
#define _RL (_RC+3*_W)

//Minized bin size to contain just one range of neighbor list
#define _NBX ((short)(_CUBEX/_RL))
#define _NBY ((short)(_CUBEY/_RL))

//-------- MDFrame Part ------------
// Set of parameters in predictor-corrector method 
#define _F02 (3./6.)
#define _F12 (251./360.)
#define _F32 (11./18.)
#define _F42 (1./6.)
#define _F52 (1./60.)      

//-------- User part ---------------
//defautl behavior of some user functions
#ifndef _USER_RESCALE
#define _USER_RESCALE() ScaleVel(0,_NP,_TR)
#endif
#ifndef _USER_INIT
#define _USER_INIT() ((void)0)
#endif
#ifndef _USER_FINISH
#define _USER_FINISH() ((void)0)
#endif
#ifndef _USER_STEP
#define _USER_STEP() ((void)0)
#endif
#ifndef _USER_SINGLE
#define _USER_SINGLE(n) ((void)0)
#endif
#ifndef _USER_PAIR
#define _USER_PAIR(i,j) ((void)0)
#endif
//=====================================================
// Debugging
//=====================================================
//kinda stupid but efficient ASSERT definitions      
#ifdef _DEBUG
#define ASSERT(p,m) if(p) ; else printf("****ASSERT FAILED:(%s,%d)",__FILE__,__LINE__), printf(m),printf("\n"),exit(1)
#define ASSERT2(p,m1,m2) if(p) ; else printf("****ASSERT FAILED:(%s,%d)",__FILE__,__LINE__ ), printf(m1,m2),printf("\n"),exit(1)
#define ASSERT3(p,m1,m2,m3) if(p) ; else printf("****ASSERT FAILED:(%s,%d)",__FILE__,__LINE__ ), printf(m1,m2,m3),printf("\n"),exit(1)
#define ASSERT4(p,m1,m2,m3,m4) if(p) ; else printf("****ASSERT FAILED:(%s,%d)",__FILE__,__LINE__ ), printf(m1,m2,m3,m4),printf("\n"),exit(1)
#else
#define ASSERT(p,m) ((void)0)
#define ASSERT2(p,m,n) ((void)0)
#define ASSERT3(p,m,n,l) ((void)0)
#define ASSERT4(p,m,n,l,k) ((void)0)
#endif
//=====================================================
// Global variables
//=====================================================
//All the global variables are small things, such as pointers or file handles

//TP is the type of particle index, if number of particles < 32767,
//  we should use `short' instead of `int'
typedef short TP;

_S double *X0,*Y0,*FX,*FY;
_S double *X1,*Y1;
// MDFrame part 
_S double *X2,*X3,*X4,*X5,*Y2,*Y3,*Y4,*Y5;
// Neib Part 
_S double *PFX,*PFY;
_S TP (*nbList)[_NNL];
_S TP *nbIndex;

_S TP (*binList)[_NBY][_NBL];
_S TP (*binIndex)[_NBY];


//---------- Multi-process control variable -------------
#define _NSEMA _NPROC+5
#define SEMA_PROP _NPROC
//share of  FIN and B1!!
#define SEMA_FIN _NPROC+1
#define SEMA_B1 _NPROC+1
#define SEMA_B2 _NPROC+2
#define SEMA_B3 _NPROC+3
#define SEMA_PASS _NPROC+4

//Start and end of nx of each process
_S short cnxs=0,cnxe=_NBX; 
_S int volatile *pstep;
_S double *pT, *pVIR, *pV;
#ifdef _RPM
_S int *fRate;
#endif
_S int ID;//Process ID, from 0 to _NPROC-1

#include "Mng.c"

//---------- Property part(extracted) -------------
#include "Property.h"
//=====================================================
// Stdlib routines
//=====================================================
#define RAND (((double)rand())/RAND_MAX)
_S double randnorm()
{
    return RAND+RAND+RAND+RAND+RAND+RAND+RAND+RAND+RAND+RAND+RAND+RAND-6;
}
//=====================================================
// Process Manage Routines
//=====================================================
_S void Block1()
{
    SemOp(SEMA_B1,-1);
    SemOp(SEMA_B1,0);
}

_S void Block3()
{
    SemOp(SEMA_B3,-1);
    SemOp(SEMA_B3,0);
}

_S void Block2()
{
    //Since everyone passed block1, it's a good time to 
    //set up final block.
    if(ID==0) SemOp(SEMA_PASS,1);
    SemOp(SEMA_B2,-1);
    SemOp(SEMA_B2,0);
}
    
_S void Finish()
{
    if(ID!=0)
    {
        //Announce finish
        SemOp(SEMA_FIN,1);
        //Get pass permittion
        SemOp(SEMA_PASS,0);
    }
    else
    {
        //Check finish for P1
        SemOp(SEMA_FIN,1-_NPROC);
        //Setup blocks
        SemOp(SEMA_B1,_NPROC);
        SemOp(SEMA_B2,_NPROC);
        SemOp(SEMA_B3,_NPROC);
        //Increase loop variable
        (*pstep)++;
        //Give permission
        SemOp(SEMA_PASS, -1);
    }
}

void ManageInit()
{
    X0=(double *)ShmAlloc(_NP,sizeof(double));
    Y0=(double *)ShmAlloc(_NP,sizeof(double));
    X1=(double *)ShmAlloc(_NP,sizeof(double));
    Y1=(double *)ShmAlloc(_NP,sizeof(double));
    X2=(double *)ShmAlloc(_NP,sizeof(double));
    Y2=(double *)ShmAlloc(_NP,sizeof(double));
    X3=(double *)ShmAlloc(_NP,sizeof(double));
    Y3=(double *)ShmAlloc(_NP,sizeof(double));
    X4=(double *)ShmAlloc(_NP,sizeof(double));
    Y4=(double *)ShmAlloc(_NP,sizeof(double));
    X5=(double *)ShmAlloc(_NP,sizeof(double));
    Y5=(double *)ShmAlloc(_NP,sizeof(double));
    PFX=(double *)ShmAlloc(_NP,sizeof(double));
    PFY=(double *)ShmAlloc(_NP,sizeof(double));
    FX=(double *)ShmAlloc(_NP,sizeof(double));
    FY=(double *)ShmAlloc(_NP,sizeof(double));
    pT=(double *)ShmAlloc(1,sizeof(double));
    pVIR=(double *)ShmAlloc(1,sizeof(double));
    pV=(double *)ShmAlloc(1,sizeof(double));
#ifdef _PROP_MESH
    meshVX=(double (*)[_MESHY])ShmAlloc(_MESHX*_MESHY,sizeof(double));
    meshVY=(double (*)[_MESHY])ShmAlloc(_MESHX*_MESHY,sizeof(double));
    mesh=(int (*)[_MESHY])ShmAlloc(_MESHX*_MESHY,sizeof(int));
#endif
#ifdef _RPM
    fRate=(int *)ShmAlloc(_NPROC,sizeof(int));
    memset(fRate,0, sizeof(int)*_NPROC);
#endif
    pstep=(int *)ShmAlloc(1,sizeof(int));
    nbList=(TP (*)[_NNL])ShmAlloc((_NP-1)*_NNL,sizeof(TP));
    nbIndex=(TP *)ShmAlloc(_NP,sizeof(TP));
    binList=(TP (*)[_NBY][_NBL])ShmAlloc(_NBX*_NBY*_NBL,sizeof(TP));
    binIndex=(TP (*)[_NBY])ShmAlloc(_NBX*_NBY,sizeof(TP));
    *pstep=0;
    *pT=*pVIR=*pV=0.0;
    SemInit(_NSEMA,_NPROC+1);
    SemOp(SEMA_B1,_NPROC);
    SemOp(SEMA_B2,_NPROC);
    SemOp(SEMA_B3,_NPROC);
}

void ManageFinish()
{
    ShmDestroy();
    SemDestroy();
}
    
//====================================================================
// BC Part
//====================================================================
//pure short hand for typing and a bit program consistency
#ifdef _2DIM
#define NearXYZ(xa,ya,za) \
{\
    rx=xa; \
    if(rx > _CUBEX/2) rx-=_CUBEX;\
    else if(rx < -_CUBEX/2) rx+=_CUBEX;\
    ry=ya; \
    if(ry > _CUBEY/2) ry-=_CUBEY;\
    else if(ry < -_CUBEY/2) ry+=_CUBEY;\
}
#else
#define NearXYZ(xa,ya,za) \
{\
    rx=xa;\
    if(rx > _CUBEX/2) rx-=_CUBEX;\
    else if(rx < -_CUBEX/2) rx+=_CUBEX;\
    ry=ya;\
    if(ry > _CUBEY/2) ry-=_CUBEY;\
    else if(ry < -_CUBEY/2) ry+=_CUBEY;\
    rz=za;\
    if(rz > _CUBEZ/2) rz-=_CUBEZ;\
    else if(rz < -_CUBEZ/2) rz+=_CUBEZ;\
}
#endif
//================================================================
// Neib part
//================================================================

//!!!!!!!!!!!This part is not optimized!!!!!!!!!!!!!!!!!!!!

//In this program, only binlist is supported now.
//Should include naive pair interaction for error checking.
#define NeibNum(n) nbIndex[n]
#define NeibVal(n,k) nbList[n][k]

#define DOCN(x,y) \
        for(l=0;l<binIndex[x][y];l++)\
        {\
            j=binList[x][y][l];\
            if(j>i)\
            {\
                NearXYZ(xi-X0[j],yi-Y0[j],zi-Z0[j]);\
                if(rx*rx+ry*ry _3(+rz*rz) < SQ(_RL))\
                {\
                    ASSERT(k<_NNL,"Out of neighbor list(Initializing)\n");\
                    nbList[i][k++]=j;\
                }\
            }\
        }

_S void NeibInit()
{
    register TP i,j,k,l;
    int n;
    short nx_m,nx_p,ny_m,ny_p;
    register double xi,yi,rx,ry _3(_ zi _ rz);
    short nx,ny _3(_ nz);
    ASSERT(_CUBEX/_NBX > _RL && _CUBEY/_NBY > _RL, "Bin too small");
    for(nx=0;nx<_NBX;nx++)
        for(ny=0;ny<_NBY;ny++)
            binIndex[nx][ny]=0;
    for(i=0;i<_NP;i++)
    {
        nx=X0[i]*_NBX/_CUBEX;
        ny=Y0[i]*_NBY/_CUBEY;
        binList[nx][ny] _3([nz])[binIndex[nx][ny]_3([nz])++]=i;
    }
    n=0;
    for(i=0;i<_NP-1;i++)
    {
        xi=PFX[i]=X0[i];
        yi=PFY[i]=Y0[i];
        _3( zi=PFZ(i)=Z0[i]);
        nx=xi*_NBX/_CUBEX;
        ny=yi*_NBY/_CUBEY;
        k=0;
        nx_m=nx-1;
        if(nx_m < 0)nx_m=_NBX-1;
        nx_p=nx+1;
        if(nx_p == _NBX) nx_p=0;
        ny_m=ny-1;
        if(ny_m < 0)ny_m=_NBY-1;
        ny_p=ny+1;
        if(ny_p == _NBY) ny_p=0;
        DOCN(nx_m,ny_m);
        DOCN(nx_m,ny);
        DOCN(nx_m,ny_p);
        DOCN(nx,ny_m);
        DOCN(nx,ny);
        DOCN(nx,ny_p);
        DOCN(nx_p,ny_m);
        DOCN(nx_p,ny);
        DOCN(nx_p,ny_p);
        nbIndex[i]=k;
        n+=k;
    }
    printf("  NEIB: Bin/DList: %dx%d bins, average %f/%d neighbors\n",
           _NBX,_NBY,((double)n)/_NP,_NNL);
}

//--------------------------------------------------------------------
// Here is a potential error, that if one particle (of small #)
// don't move for a long time, it's neighbor list will eventually
// occupied by the empty slots(-1), though, this will increase
// 12%?? percent of performance of the program.
#define CKNB(x,y,lb) \
    for(l=0;l<binIndex[x][y];l++)\
    {\
        j=binList[x][y][l];\
        if(j<i)\
        {\
             /*=====UpList*/ \
            NearXYZ(xi-X0[j],yi-Y0[j],zi-Z0[j]);\
            r2=rx*rx+ry*ry _3(+rz*rz);\
            if(r2 > _RL*_RL)\
            {\
                if(r2 <= (_RL+3*_W)*(_RL+3*_W))\
                    for(k=0;k<nbIndex[j];k++)\
                        if(nbList[j][k] == i)\
                            nbList[j][k]=-1;\
            }\
            else if(r2 >= (_RL-3*_W)*(_RL-3*_W))\
            {\
                for(k=0;k<nbIndex[j];k++)\
                    if(nbList[j][k] == i) goto lb;\
                ASSERT(k<_NNL,"Out of neighbor list(Running)\n");\
                nbList[j][k++]=i;\
                nbIndex[j]=k;\
            }\
        }\
        else if(j>i)\
        {\
             /*=====SelfList*/\
            NearXYZ(xi-X0[j],yi-Y0[j],zi-Z0[j]);\
            if(rx*rx+ry*ry _3(+rz*rz) < SQ(_RL))\
                nbList[i][ksum++]=j;\
        }\
    lb:;\
    }

_S void NeibBinFlashAtom(TP i)
{
    register TP j,k,l,ksum;
    register short nx,ny;
    short nx_m,nx_p,ny_m,ny_p;
    register double xi,yi,rx,ry _3(_ zi _ rz), r2;
    xi=PFX[i]=X0[i];
    yi=PFY[i]=Y0[i];
    _3( zi=PFZ(i)=Z0[i]);
    nx=xi*_NBX/_CUBEX;
    ny=yi*_NBY/_CUBEY;
    ksum=0;
    nx_m=nx-1;
    if(nx_m < 0)nx_m=_NBX-1;
    nx_p=nx+1;
    if(nx_p == _NBX) nx_p=0;
    ny_m=ny-1;
    if(ny_m < 0)ny_m=_NBY-1;
    ny_p=ny+1;
    if(ny_p == _NBY) ny_p=0;
    CKNB(nx_m,ny_m,l1);
    CKNB(nx_m,ny,l2);
    CKNB(nx_m,ny_p,l3);
    CKNB(nx,ny_m,l4);
    CKNB(nx,ny,l5);
    CKNB(nx,ny_p,l6);
    CKNB(nx_p,ny_m,l7);
    CKNB(nx_p,ny,l8);
    CKNB(nx_p,ny_p,l9);
    nbIndex[i]=ksum;
    //printf("Flashing %d, neib %d\n",i,ksum);
}

_S void NeibStep()
{
    register TP i,j,l;
    register short nx,ny _3(_ nz),cx,cy;
    for(nx=cnxs;nx<cnxe;nx++)
    {
        if(nx==cnxs) SemLock(ID);
        else if(nx==cnxe-1) SemLock((ID+1)%_NPROC);
        
        for(ny=0;ny<_NBY;ny++) for(l=0;l<binIndex[nx][ny];l++) 
            //for(i=0;i<_NP;i++)
        {
            i=binList[nx][ny][l];
            cx=X0[i]*_NBX/_CUBEX;
            cy=Y0[i]*_NBY/_CUBEY;
            _3(cz=Z0[i]*_NBZ/_CUBEZ);
            ASSERT(cx<_NBX && cx>=0 && cy<_NBY && cy >=0,
                   "NeibStep, out of frame");
            if(nx!=cx || ny!=cy)
            {
#ifdef _RPM 
                //There is a flip here, reflecting a particle only change it's
                //odd order derivative of position, so it will pass the RPM
                //plane and pass back immediately.
                if(cx==nx-1 && nx>=_RPMIN && nx<=_RPMAX)
                {
                    if(nx==_RPMIN) fRate[ID]--;
                    ASSERT(X1[i] < 0,"How can you pass with negative speed?");
                    if(RAND <= _RPM)
                    {
                        X1[i]=-X1[i];
                        X3[i]=-X3[i];
                        X5[i]=-X5[i];
                    }
                }
                else if(cx==nx+1 && nx==_RPMIN) fRate[ID]++;
#endif
                //$memcpy
                for(j=l;j<binIndex[nx][ny];j++)
                    binList[nx][ny][j]=binList[nx][ny][j+1];
                binIndex[nx][ny]--;
                binList[cx][cy]_3([cz])[binIndex[cx][cy]++]=i;
                ASSERT(binIndex[cx][cy] < _NBL, "Out of Binlist");
            }
            if(SQ(PFX[i]-X0[i])+SQ(PFY[i]-Y0[i])_3(+SQ(PFZ[i]-Z0[i]))>=SQ(_W))
                    NeibBinFlashAtom(i);
        }
        if(nx==cnxs) SemUnlock(ID);
        else if(nx==cnxe-1) SemUnlock((ID+1)%_NPROC);
    }        
}
//--------------------------------------------------------------------

_S void NeibFinish()
{
}
//================================================================
// MDFrame Part
//================================================================
_S void MPInit()
{
    printf("  MDFrame: Using Gear's 5th order pred/corr method\n");
}

//-------   The Main subroutine of MDFrame.F -------------------------
_S void MoveParticle()
{
    register TP i,l;
    short nx,ny;
    //----------- Predictor
    for(nx=cnxs;nx<cnxe;nx++) for(ny=0;ny<_NBY;ny++)
        for(l=0;l<binIndex[nx][ny];l++) 
    //for(i=0;i<_NP;i++)
    {
        i=binList[nx][ny][l];
        X0[i]+=X1[i]+X2[i]+X3[i]+X4[i]+X5[i];
        X1[i]+=2*X2[i]+3*X3[i]+4*X4[i]+5*X5[i];
        X2[i]+=3*X3[i]+6*X4[i]+10*X5[i];
        X3[i]+=4*X4[i]+10*X5[i];
        X4[i]+=5*X5[i];
        FX[i]=0;

        Y0[i]+=Y1[i]+Y2[i]+Y3[i]+Y4[i]+Y5[i];
        Y1[i]+=2*Y2[i]+3*Y3[i]+4*Y4[i]+5*Y5[i];
        Y2[i]+=3*Y3[i]+6*Y4[i]+10*Y5[i];
        Y3[i]+=4*Y4[i]+10*Y5[i];
        Y4[i]+=5*Y5[i];
        FY[i]=0;
        
        _3( Z0[i]+=Z1[i]+Z2[i]+Z3[i]+Z4[i]+Z5[i];
            Z1[i]+=2*Z2[i]+3*Z3[i]+4*Z4[i]+5*Z5[i];
            Z2[i]+=3*Z3[i]+6*Z4[i]+10*Z5[i];
            Z3[i]+=4*Z4[i]+10*Z5[i];
            Z4[i]+=5*Z5[i];
            FZ[i]=0);
    }
    Block3();
    //------------ Interactions
    for(nx=cnxs;nx<cnxe;nx++)
    {
        if(nx==cnxs) SemLock(ID);
        else if(nx==cnxe-1) SemLock((ID+1)%_NPROC);
        
        for(ny=0;ny<_NBY;ny++) for(l=0;l<binIndex[nx][ny];l++)
            //for(i=0;i<_NP-1;i++)
        {
            register TP j,k;
            register double r2,rx,ry,xi,yi,fxi,fyi _3(_ zi _ fzi _ rz);
            i=binList[nx][ny][l];
            fxi=fyi= _3(fzi=) 0;
            xi=X0[i];
            yi=Y0[i];
            _3( zi=Z0[i]);
            for(k=0;k<NeibNum(i);k++)
            {
                j=NeibVal(i,k);
                if(j != -1)
                {
                    //#BC entry
                    NearXYZ( xi-X0[j],yi-Y0[j],zi-Z0[j]);
                    r2=rx*rx+ry*ry _3(+rz*rz);
                    if(r2 <= _RC*_RC)
                    {
                        //#Pot Entry
                        register double r6,ff;
                        ASSERT4(r2 > 0.65, "Too Near:(%d %d)%f",i,j,r2);
                        r6=r2*r2*r2;
                        ff=24*(2/r6-1)/(r2*r6);
                        fxi += ff*rx;
                        FX[j] -= ff*rx;
                        fyi += ff*ry;
                        FY[j] -= ff*ry;
                        _3( fzi += ff*rz;
                            FZ[j] -= ff*rz);
                        //#Prop Entry
                        if(*pstep%_NDUMP == 0)
                        {
                            V+=4/r6*(1/r6-1)-4*(1/(_RC*_RC*_RC*_RC*_RC*_RC*
                                                   _RC*_RC*_RC*_RC*_RC*_RC)
                                                -1/(_RC*_RC*_RC*_RC*_RC*_RC));
                            VIR+=ff*r2;
                        }
                        _USER_PAIR(i,j);
                    }
                }
            }
//#if _POT_SINGLE
#ifdef _CR
            //Center Region
            if(xi < _CR_X1+_RL && xi > _CR_X0-_RL &&
               yi < _CR_Y1+_RL && yi > _CR_Y0-_RL)
            {
                for(j=_CR_X0;j<=_CR_X1;j++) for(k=_CR_Y0;k<=_CR_Y1;k++)
                {
                    rx=xi-j;
                    ry=yi-k;
                    r2=rx*rx+ry*ry;
                    if(r2 <= _RC*_RC)
                    {
                        //#Pot Entry
                        register double r6,ff;
                        ASSERT4(r2 > 0.65, "Too Near:(%d %d)%f",i,j,r2);
                        r6=r2*r2*r2;
                        ff=24*(2/r6-1)/(r2*r6);
                        fxi += ff*rx;
                        fyi += ff*ry;
                        //#Prop Entry
                        if(*pstep%_NDUMP == 0)
                        {
                            V+=4/r6*(1/r6-1)-4*(1/(_RC*_RC*_RC*_RC*_RC*_RC*
                                                   _RC*_RC*_RC*_RC*_RC*_RC)
                                                -1/(_RC*_RC*_RC*_RC*_RC*_RC));
                            VIR+=ff*r2;
                        }
                    }
                }
            }
#endif
//#endif
            FX[i]+=fxi;
            FY[i]+=fyi;
            _3( FZ[i]+=fzi);
        }
        if(nx==cnxs) SemUnlock(ID);
        else if(nx==cnxe-1) SemUnlock((ID+1)%_NPROC);
    }
    Block2();
    //-------------- Corrector
    for(nx=cnxs;nx<cnxe;nx++) for(ny=0;ny<_NBY;ny++)
        for(l=0;l<binIndex[nx][ny];l++) 
    //for(i=0;i<_NP;i++)
    {
        register double xi,yi _3(_ zi);
        i=binList[nx][ny][l];
        xi=X2[i]-FX[i]*(_DT*_DT/2);
        X0[i]-=xi*_F02;
        if(X0[i] >= _CUBEX) X0[i]-=_CUBEX;
        else if(X0[i] < 0) X0[i]+=_CUBEX;
        X1[i]-=xi*_F12;
        X2[i]-=xi;
        X3[i]-=xi*_F32;
        X4[i]-=xi*_F42;
        X5[i]-=xi*_F52;
        
        yi=Y2[i]-FY[i]*(_DT*_DT/2);
        Y0[i]-=yi*_F02;
        if(Y0[i] >= _CUBEY) Y0[i]-=_CUBEY;
        else if(Y0[i] < 0) Y0[i]+=_CUBEY;
        Y1[i]-=yi*_F12;
        Y2[i]-=yi;
        Y3[i]-=yi*_F32;
        Y4[i]-=yi*_F42;
        Y5[i]-=yi*_F52;
        
        _3( zi=Z2[i]-FZ(i)*(_DT*_DT/2);
            Z0[i]-=zi*_F02;
            if(Z0[i] >= _CUBEZ) Z0[i]-=_CUBEZ;
            else if(Z0[i] < 0) Z0[i]+=_CUBEZ;
            Z1[i]-=zi*_F12;
            Z2[i]-=zi;
            Z3[i]-=zi*_F32;
            Z4[i]-=zi*_F42;
            Z5[i]-=zi*_F52);
            
        ASSERT3(X0[i] < _CUBEX, "X out:%d,%f",i,X0[i]);
        ASSERT3(Y0[i] < _CUBEY, "Y out:%d,%f",i,Y0[i]);
        _3(ASSERT3(Z0[i] < _CUBEZ, "Z out:%d,%f",i,Z0[i]));
        //#Prop Entry
        if(*pstep%_NDUMP == 0)
        {
            T+=X1[i]*X1[i]+Y1[i]*Y1[i] _3( +Z1[i]*Z1[i]);
        }
        PropSingle(i);
        _USER_SINGLE(i);
    }
}

_S void MPFinish()
{
}
//=====================================================
// Frame routines
//=====================================================
//------------- Initialization helper starts here---------------------
//Used for place solid in good(viewable) position
//#define D_SC 0 
//#define D_FCC 0
//#define D_BCC 0

#ifndef _2DIM
//---------- InitSC: Init from particle# bn to en as SC---------------
_S void InitSC(double bx,double by,double bz,
               double sx,double sy,double sz,TP bn,TP sn)
{
    TP n;
    short nx,ny,nz;
    double ddx,ddy,ddz,rx,ry,rz;
    
    nx=cbrt(sn*sx*sx/sy/sz)+1;
    ny=cbrt(sn*sy*sy/sz/sx);
    nz=cbrt(sn*sz*sz/sx/sy);
    if(nx*ny*nz < sn)
    {
        ny++;
        if(nx*ny*nz < sn) nz++;
    }
    ddx=(sx-1e-9)/nx;
    ddy=(sy-1e-9)/ny;
    ddz=(sz-1e-9)/nz;

    n=bn;
    for(rz=bz;rz<=sz+bz-ddz;rz+=ddz)
        for(ry=by;ry<=sy+by-ddy;ry+=ddy)
            for(rx=bx;rx<=sx+bx-ddx;rx+=ddx)
            {
                X0[n]=rx;
                Y0[n]=ry;
                Z0[n]=rz;
                n++;
                if(n >= sn+bn) return;
            }
    printf("***Error: Frame: Not enough space to init SC\n");
}
//--------------------------------------------------------------------
_S void InitFCC(double bx,double by,double bz,
                double sx,double sy,double sz,TP bn,TP sn)
{
    TP n;
    short nx,ny,nz;
    double ddx,ddy,ddz,rx,ry,rz;
    
    nx=cbrt(sn/4*sx*sx/sy/sz)+1;
    ny=cbrt(sn/4*sy*sy/sz/sx);
    nz=cbrt(sn/4*sz*sz/sx/sy);
    if(nx*ny*nz*4 < sn)
    {
        ny++;
        if(nx*ny*nz*4 < sn)nz=nz+1;
    }
    ddx=(sx-1e-9)/nx;
    ddy=(sy-1e-9)/ny;
    ddz=(sz-1e-9)/nz;
    
    n=bn;
    rz=0;
    for(rz=0;rz<=sz-ddz;rz+=ddz)
        for(ry=0;ry<=sy-ddy;ry+=ddy)
            for(rx=0;rx<=sx-ddx;rx+=ddx)
            {
                X0[n]=rx+bx+ddx/4;
                Y0[n]=ry+by+ddy/4;
                Z0[n]=rz+bz+ddz/4;
                n++;
                if(n >= sn+bn) return;
                X0[n]=rx+bx+ddx/4;
                Y0[n]=ry+by+ddy/4+ddy/2;
                Z0[n]=rz+bz+ddz/4+ddz/2;
                n++;
                if(n >= sn+bn) return;
                X0[n]=rx+bx+ddx/4+ddx/2;
                Y0[n]=ry+by+ddy/4;
                Z0[n]=rz+bz+ddz/4+ddz/2;
                n++;
                if(n >= sn+bn) return;
                X0[n]=rx+bx+ddx/4+ddx/2;
                Y0[n]=ry+by+ddy/4+ddy/2;
                Z0[n]=rz+bz+ddz/4;
                n++;
                if(n >= sn+bn) return;
            }
    printf("***Error: Frame: Not enough space to init FCC\n");
}
#else
//----------  InitSC2D: Init from particle# bn to en as SC for 2D-----
//This initialization is not good for _CR defined, only temperarily used
_S void InitSC2D(double bx,double by,double sx,double sy,TP bn,TP sn)
{
    TP n;
    short nx,ny;
    double ddx,ddy,rx,ry;
#ifdef _CR
    //$TMP
    sn+=400;
    //END$TMP
#endif
    ddx=sqrt(sn*sx/sy);
    nx=(ddx+1);
    ny=(sn/ddx);
    if(nx*ny < sn) ny++;
    ddx=(sx-1e-9)/nx;
    ddy=(sy-1e-9)/ny;
#ifdef _CR
    //$TMP
    sn-=400;
    nx+=4;
    ny+=4;
    //END$TMP
#endif
    n=bn;
    for(ry=0;ry<=sy-ddy;ry+=ddy)
        for(rx=0;rx<=sx-ddx;rx+=ddx)
        {
            X0[n]=rx+bx+ddx/2;
            Y0[n]=ry+by+ddy/2;
#ifdef _CR
            //$TMP
            if(X0[n] <= _CR_X1+1 && X0[n] >= _CR_X0-1 &&
               Y0[n] <= _CR_Y1+1 && Y0[n] >= _CR_Y0-1)
                continue;
            //END$TMP
#endif
            n++;
            if(n >= sn+bn) return;
        }
    printf("***Error: Frame: Not enough space to init SC\n");
}
#endif
//---------- Initialize velocity according to the temperature---------
_S void InitVel(TP bn, TP en, double tp)
{
    TP i;
    double s,px,py _3(_  pz);
    s=sqrt(tp);
    px=0;
    py=0;
    _3(pz=0);
    for(i=bn;i<en;i++)
    {
        X1[i]=randnorm()*s*_DT;
        px+=X1[i];
        Y1[i]=randnorm()*s*_DT;
        py+=Y1[i];
        _3( Z1[i]=randnorm()*s*_DT);
        _3( pz+=Z1[i]);
    }
    px/=(en-bn);
    py/=(en-bn);
    _3( pz/=(en-bn));
    for(i=bn;i<en;i++)
    {
        X1[i]-=px;
        Y1[i]-=py;
        _3( Z1[i]-=pz);
    }
}
//------------ Rescale velocity according to the temperature--------------
_S void ScaleVel(TP bn,TP en,double tp)
{
    double t;
    TP i;
    t=0;
    for(i=bn;i<en;i++)
        t+=X1[i]*X1[i]+Y1[i]*Y1[i] _3(  +Z1[i]*Z1[i]);
    t=sqrt(tp*_DIMINT*(en-bn)/t)*_DT;
    for(i=bn;i<en;i++)
    {
        X1[i]*=t;
        Y1[i]*=t;
        _3( Z1[i]*=t);
    }
}
#ifndef _2DIM      
#define DoInit() (InitFCC(0,0,0,_CUBEX,_CUBEY,_CUBEZ,0,_NP),InitVel(0,_NP,_TR))
#else
#define DoInit() (InitSC2D(0,0,_CUBEX,_CUBEY,0,_NP),InitVel(0,_NP,_TR))
#endif
//---------- The main running subroutine ---------------------
_S void Run()
{
    int iTMP,handle;
#ifdef _PROP_PARTICLE
    FILE *part_fp;
#endif
    register int i;
    double dTMP;
    ManageInit();
    printf("======================================================\n"
    //---- Dumping      
           "Atomic Simulation Library(C/C++) (Version %s) : %s\n",
           _VERSION, _APPLICATION);
    printf(" Initializing...\n");
    //---- Init
#ifndef _LCFG
    DoInit();
#else 
    handle=open(_LCFG_READ,O_RDONLY,0600);
    if(handle < 0)
    {
        printf("****Warning:Can't open %s to read CFG, continue w/o CFG\n",
               _LCFG_READ);
        DoInit();
    }
    else
    {
        read(handle,&iTMP,sizeof(int));
        if(iTMP!=_DIMINT) printf("Warning: DIMINT not same\n");
        read(handle,&iTMP,sizeof(int));
        if(iTMP!=_NP) printf("Warning: NP not same\n");
        read(handle,&dTMP,sizeof(double));
        if(dTMP!=_TR) printf("Warning: TR not same\n");
        read(handle,&dTMP,sizeof(double));
        if(dTMP!=_CUBEX) printf("Warning: CUBEX not same\n");
        read(handle,&dTMP,sizeof(double));
        if(dTMP!=_CUBEY) printf("Warning: CUBEY not same\n");
        _3( read(handle,&dTMP,sizeof(double));
            if(dTMP!=_CUBEZ) printf("Warning: CUBEZ not same\n") );
        read(handle,X0,_NP*sizeof(double));
        read(handle,Y0,_NP*sizeof(double));
        _3( read(handle,&Z0,_NP*sizeof(double)) );
        read(handle,X1,_NP*sizeof(double));
        read(handle,Y1,_NP*sizeof(double));
        _3( read(handle,&Z1,_NP*sizeof(double)) );
        for(i=0;i<_NP;i++)
        {
            X1[i]*=_DT;
            Y1[i]*=_DT;
            _3( Z1[i]*=_DT );
        }
        close(handle);
        printf("  CFG: Loaded from %s\n",_LCFG_READ);
    }
#endif 
    printf(" Basic parameters: particle number : %d\n"
           "    Density : %f, Temperature: %f\n"
           "    Cube : %f X %f"_3(" X %f")"\n"
           ,_NP,(double)_DR,(double)_TR,
           (double)_CUBEX, (double)_CUBEY _3(_ (double)_CUBEZ));

    NeibInit();
    _USER_INIT();//Notice the position of UserInit()
    MPInit();
    PropInit();
    //==== Spawning children
    cnxe=0;
#if _NPROC == 1
    printf(" One process serial code\n");
#else
    printf(" Spawning %d children...\n",_NPROC);
    iTMP=_NBX%_NPROC;
    for(ID=1;ID<_NPROC;ID++)
    {
        cnxs=cnxe;
        cnxe=cnxs+(short)((_NBX)/_NPROC);//$
        if(iTMP > 0)
        {
            cnxe++;
            iTMP--;
        }
        if(fork()==0) goto l_Dothings;
        printf("    Child %d:BinX(%d--%d)\n",ID,cnxs,cnxe-1);
    }
#endif
    cnxs=cnxe;
    cnxe=_NBX;
    ID=0;
    printf("    Master:BinX(%d--%d)\n",cnxs,cnxe-1);
    //==== Start cycling
    printf(" Start running...\n");
#ifdef _PROP_PARTICLE
    part_fp=fopen(_PROP_PARTICLE,"wt");
    if(part_fp == 0)
        printf("****Warning:Can't open %s to write particle data\n",
               _PROP_PARTICLE);
#endif
 l_Dothings:
#if _NRESCALE
    //---- Rescale
    if(ID==0) printf("  Start rescaling ...\n");
    for(*pstep=0;*pstep<_NRESCALE;)
    {
        NeibStep();
        Block1();
        if(ID==0) _USER_RESCALE();
        MoveParticle();
        PropStep();
        _USER_STEP();
        Finish();
        if(*pstep % _NDUMP==0) if(ID==0)
        {
            printf("%d\t%f\t%f\t%.8f\n",*pstep, *pT/_NP/_DT/_DT/_DIMINT,
                   (*pVIR+*pT/_DT/_DT)/_DIMINT/_V, (*pT/_DT/_DT/2+*pV)/_NP);
            *pT=*pVIR=*pV=0;
        }
    }
    if(ID==0) printf("  End rescaling.\n");
#endif */
//---- Normal Run
    //for(step=_NRESCALE;step<_NSTEP;step++)
    for(*pstep=_NRESCALE;*pstep<_NSTEP;)
    {
        NeibStep();
        Block1();
        MoveParticle();
        PropStep();
        _USER_STEP();
        Finish();
        if(ID==0)
        {
#ifdef _PROP_PARTICLE
            if(*pstep % _NCORR==0)
                for(i=0;i<_NP;i++) fprintf(part_fp,"%f %f %f %f\n",
                                           X0[i],Y0[i],X1[i]/_DT,Y1[i]/_DT);
#endif
            if(*pstep % _NDUMP==0) 
            {
                printf("%d\t%f\t%f\t%.8f\n",*pstep, *pT/_NP/_DT/_DT/_DIMINT,
                       (*pVIR+*pT/_DT/_DT)/_DIMINT/_V,
                       (*pT/_DT/_DT/2+*pV)/_NP);
                *pT=*pVIR=*pV=0;
            }
        }
    }
    if(ID==0)
    {
        //==== Finished      
        printf(" Finished! Exiting...\n");
#ifdef _PROP_PARTICLE
        fclose(part_fp);
#endif
#ifdef _RPM
        //print the flow rate
        for(i=0,iTMP=0;i<_NPROC;i++)
            iTMP+=fRate[i];
        printf("  Flow rate:%d particles/%d timesteps(%f/reduced time)\n",
               iTMP,_NSTEP,iTMP/_DT/_NSTEP);
#endif
        handle=open(_LCFG_WRITE,O_WRONLY|O_CREAT,0600);
        if(handle < 0) printf("****Error:Can't open %s to write CFG\n",
                              _LCFG_WRITE),exit(1);
        iTMP=_DIMINT;
        write(handle,&iTMP,sizeof(int));
        iTMP=_NP;
        write(handle,&iTMP,sizeof(int));
        dTMP=_TR;
        write(handle,&dTMP,sizeof(double));
        dTMP=_CUBEX;
        write(handle,&dTMP,sizeof(double));
        dTMP=_CUBEY;
        write(handle,&dTMP,sizeof(double));
        _3( dTMP=_CUBEZ;
            write(handle,&dTMP,sizeof(double)));
        write(handle,X0,_NP*sizeof(double));
        write(handle,Y0,_NP*sizeof(double));
        _3( write(handle,Z0,_NP*sizeof(double)) );
        for(i=0;i<_NP;i++)
        {
            X1[i]/=_DT;
            Y1[i]/=_DT;
            _3( Z1[i]/=_DT );
        }
        write(handle,X1,_NP*sizeof(double));
        write(handle,Y1,_NP*sizeof(double));
        _3( write(handle,Z1,_NP*sizeof(double)) );
        _USER_FINISH();
        PropFinish();
        MPFinish();
        NeibFinish();
        close(handle);
        printf("  CFG: Wrote to %s\n",_LCFG_WRITE);
        printf("======================================================\n");
        ManageFinish();
    }
}

//=====================================================
// Main
//=====================================================
#ifndef _USER_MAIN
int main()
{
    Run();
    return 0;
}
#endif