Forces.h

#ifndef FORCES_H
#define FORCES_H
#include <VarData.h>
#include <Cubo.h>
#include <SingProc.h>
#include <time.h>
#ifdef USE_GL
#include <GL/glut.h>
#endif//USE_GL

enum SYS_SHAPE{
  SYS_1D =     0x0001,
  SYS_2D =     0x0002,
  SYS_3D =     0x0004,
  SYS_LEAVES = 0x0008,
  SYS_RIGID  = 0x0010,
  SYS_STALK  = 0x0020,
  SYS_MD     = 0x0040,
  SYS_MC     = 0x0080,
  SYS_WIDOM  = 0x0100,
  SYS_TRIAL  = 0x0200,
  SYS_PORE   = 0x0400,
  SYS_ROD   = 0x0800,
  SYS_ELECTRO   = 0x1000,
};
enum CALC_MODE{
  CALC_NVT     = 0x00001,
  CALC_mVT     = 0x00002,
  CALC_SOLVE   = 0x00004,
  CALC_PAIR    = 0x00008,
  CALC_HARM    = 0x00010,
  CALC_LJ39    = 0x00020,
  CALC_LJ      = 0x00040,
  CALC_DENS_CH = 0x00080,
  CALC_NcVT    = 0x00100,
  CALC_mcVT    = 0x00200,
  CALC_3d      = 0x00400,
  CALC_2d      = 0x00800,
  CALC_DENS    = 0x01000,
  CALC_CONF_BIAS = 0x02000,
  CALC_BIL_BIAS  = 0x04000,
  CALC_SPH_BIAS  = 0x08000,
  CALC_STEP  = 0x10000,
  CALC_ELECTRO  = 0x20000,
};
enum THERM_MODE{
  THERM_NO      = 0x0001,
  THERM_LANG    = 0x0002,
  THERM_AND     = 0x0003,
  THERM_BERE    = 0x0004,
};
enum FIT_TYPE{
  FIT_SPLINE3 = 0,
  FIT_SPLINE4 = 1,
  FIT_PARAB2  = 2,
  FIT_CUBIC   = 3,
  FIT_DERIV   = -1,
  FIT_FORTH   = 4,
  FIT_POLY    = 5,
  FIT_BSPLINE = 6,
};
enum IntMethod{
  INT_MD = 0,
  INT_DIFF = 1,
  INT_MC = 2,
};
enum ALLOCATED{
  ALL_FORCES = 0x001,
  ALL_METR   = 0x002,
  ALL_POT    = 0x004,
  ALL_MC     = 0x008,
  ALL_MD     = 0x010,
  ALL_SPLINE = 0x020,
  ALL_CYL    = 0x040,
  ALL_FFIELD = 0x080,
  ALL_TENS   = 0x100,
  ALL_DENS   = 0x200,
  ALL_BIAS   = 0x400,
  ALL_MATRIX = 0x800,
};
typedef struct{
  double Lap;
  double SLap;
  double El[3];
  double Ext;
  double LJ;
  double Cont;
  double Elong[3];
  double LJMin;
  double LJCutOff;
  double CutOff2;
  double BaseLine;
  double PotThr;
  double ForThr;
  double DistThr;
} KFORCES;
typedef struct{
  double Dir[3];
  /* /// Helfrich */
  /* double Hel[3]; */
  double Ext[3];
  /* /// Elastic */
  /* double El[3]; */
  /* /// Lennard-Jones */
  /* double LJ[3]; */
} FORCES;
typedef struct{
  int NComp;
  int NDim;
  int NSlab;
  int CalcMode;
  double **Pre;
  double **Dens;
  double PreTot[3];
  double RefPos[3];
  double Edge[3];
  int Wrap[3];
  double EdgeInv[3];
}TENS;
//typedef DdDoubleLoop DomDec;
typedef DdLinkedList DomDec;
//typedef DdFixedSize DomDec;
//typedef DdArray DomDec;
class Forces : public VarData{
 private:
  char ConfFile[60];
  time_t InitTime;
  time_t CurrTime;
  double *Dens2;
  double *Dens3;
  double *PTab;
  double *FTab;
  double *MTab;
  double OldNrgSys;
  double *OldNrgBead;
  double *OldNrgCh;
  double **OldPos;
  double *LocDens2;
  double *LocDens3;
  double *FirstBeadDistr;
  double **BondPosBias;
  double *CumProbBias;
  double BorderBias[2];
  int NBin;
  int NGrid;
  int NTrialBias;
  int BoundCond[6];
  int PeriodicImage[3];
  FILE *StatFile1;
  FILE *StatFile2 ;
  typedef double(Forces::*CALC_NRG)(int p,double *Pot);
  CALC_NRG NrgBead;
  CALC_NRG NrgCh;
  //static CALC_NRG ChooseCalcMode(int Mode);
  void ChooseCalcMode(int Mode);
  typedef double(Forces::*CALC_POT)(double Dist,int t1,int t2,double *Pot);
  CALC_POT CalcPot;
  //static CALC_POT ChoosePot(int Mode);
  void ChoosePot(int Mode);
  int NTab;
  int DynFlag;
  typedef double(Forces::*STEP_CONF)();
  typedef double(Forces::*STEP_SIM)();
  STEP_CONF ConfigSys();
  STEP_SIM CalcUpdate();
  Matrice *IntMatrix;
  public:
  //-----------Forces.cpp---------------------
  Forces(int argc,char **argv,int NPart,char *ConfFile);
  Forces(int argc,char **argv,char *ConfFileExt,char *Snapshot);
  void Shout(const char * s, ...);
  void AllocMethod();
  void PrepareSys();
  void PrepareParallel(int argc,char **argv);
  ~Forces();
  void Info();
  int ReadConfDinamica(char *File);
  void InitConst();
  int ReSetNPart(int NewNPart);
  int ReSetNChain(int NewNChain);
  void ReSetNPCh(int NewNPCh);
  void ReOpen(char *FName,int Bf);
  void FillMatrix();
  //-----------ForcesIntegration.cpp----------
  void StudySys();
  int MinHelfrich();
  int DynIntegration();
  void Dynamics();
  void Solve();
  void SolveLinks();
  void SolveLinksSparse();
  void SolveLinksIterative();
  void SolveRod();
  void SolveLeaves();
  int Update();
  void VelVerletRigid();
  void VelVerletRigid2();
  int IfMetropolis(double ExpArg,double Weight);
  int InsertBead(int p);
  void IgnoreCh(int c);
  void RemChFromSys(int c);
  void SaveCh(int c);
  void ReInsertCh(int c);
  void ConsiderCh(int c);
  double InsertCh(int c);
  double InsertRest(int pCurr,int StartPos);
  double RemoveChBias(int c);
  double InsertChBias(int c);
  double WeightSetBond(int p,int t);
  double CreateSetBond(int p,int t);
  int MoveBead(int p);
  int TryInsert();
  int TryRemove();
  int TryMove();
  int TryMoveCh();
  int TryInsertCh();
  int TryRemoveCh();
  int TryInsertChBias();
  int TryRemoveChBias();
  void WidomInsert(double *NrgDiff);
  void WidomRemove(double *NrgDiff,int p);
  void WidomInsertCh(double *NrgDiff);
  void WidomRemoveCh(double *NrgDiff,int c);
  void WidomBiasChIn(double *Weight);
  void WidomBiasChOut(double *Weight,int c);
  void VelVerlet1();
  void LangevinTherm();
  void AndersenTherm();
  void BerendsenTherm();
  void NoTherm(){};
  void VelVerlet2();
  //static CALC_NRG ChooseCalcMode(int Mode);
  void ChooseThermostat(int Mode);
  typedef void (Forces::*CALC_THERM)();
  CALC_THERM CalcTherm;
  void ApplyTherm(){return (*this.*CalcTherm)();};
  //-----------ForcesForceField.cpp-----------
  int ForceFieldLine();
  int SumSomeForces(int HowMany);
  int ForceFieldLeaves();
  int ForceFieldBulk();
  void ForceFieldRod();
  void ForceFieldRigid();
  void CalcForcesDensFunc();
  void PrintForce();
  void GetForceField();
  void TabForceAlloc(int NTabExt);
  void TabPot();
  void NanoInteraction();
  double NanoNrg(int p);
  double NanoNrg(double *Pos,int t);
  void DefForceParam();
  void DefNanoForceParam();
  double RigidLJ(int nNano,double Dist,double *Pot,double Sign);
  double RigidHamaker(int n,double Dist,double *Pot,double Sign);
  double RigidCoulomb(int nNano,double Dist,double *Pot,double Sign);
  double RigidDistanceRad(int n,int nn,double *dr);
  double RigidDistanceAxis(int n,int nn,double *dr);
  // Point to the potential 
  void PointShape(int iShape);
  /* /// Data type for distance/field functions */
  /* typedef double(Forces::*NANO_FORCE)(double Dist2,int n,int p); */
  /* /// Pointer to a distance/field function */
  /* NANO_FORCE Nano_Force; */
  /* /// Pointer to a generic function */
  /* double NanoForce(double Dist2,int n,int p){return (*this.*Nano_Force)(Dist2,n,p);} */
  double Harmonic(double Dist2,int t1,int t2,double *Pot);
  double StepPot(double Dist2,int t1,int t2,double *Pot);
  double ElectroPot(double Dist2,int t1,int t2,double *Pot);
  double LJ39(double Dist2,int t1,int t2,double *Pot);
  double LJPot(double Dist2,int t1,int t2,double *Pot);
  double Potential(double Dist,int t1,int t2,double *Pot){return (*this.*CalcPot)(Dist,t1,t2,Pot);};  
  //-----------ForcesCalcNrg.cpp-----------
  double CalcTotNrgCh();
  double CalcTotNrgBead();
  double CalcNrgBead(int p,double *Pot){return (*this.*NrgBead)(p,Pot);};
  double CalcNrgCh(int c,double *Pot){return (*this.*NrgCh)(c,Pot);};
  double NrgChBondDens(int c,double *Pot);
  double NrgChDens(int c,double *Pot);
  double CalcPairwise(int p,double *Pot);
  double CalcPairwiseCh(int c,double *Pot);
  double CalcBending(int p);
  double CalcSpring(int p);
  double CalcBonded(int p,double *Pot);
  double CalcBendingGhost(double *Pos,int pExt);
  double CalcBondedCh(int c,double *Pot);
  /* ///Calculate the spring, bending and non bonded interactions */
  /* double CalcNrgCh(int c,double *Pot); */
  void CalcNrgBeadDensFunc();
  double CalcNrgBeadDensFunc(int p,double *Pot);
  double CheckDomDec(int p);
  void CheckPairList();
  double SumForcesMD();
  double Wei3(const double r, const double a);
  double DerWei3(const double r, const double a);
  double Wei2(const double r, const double b);
  double DerWei2(const double r, const double b);
  double DensFuncNrgGhost(double *Pos,int p1,int t1);
  double DensFuncNrgGhostInternal(double *Pos,int p1,int t1);
  double DensFuncNrgBead(int p1);
  double DensFuncNrgCh(int c,double *Pot);
  double DensFuncNrgChAv(int c);
  double DensFuncNrgChInternal(int c);
  double DensFuncNrgSys();
  void CalcDens(int pInit,int pEnd);
  void ClearDens();
  double NrgStep(int p);
  double NrgStepCh(int c,double *Pot);
  double NrgElectro(int p);
  int AddDens(int pInit,int pEnd);
  int RemDens(int pInit,int pEnd);
  double SumDens(int pInit,int pEnd);
  int ListNeiCell(int p,double *Pos,int *NeiList){
    return Pc->GetNei(Pos,NeiList);
    //return Pc->GetCell(p,NeiList);
  };
  //-----------ForcesBoundary.cpp-------------
  void PullBead();
  void PushBead();
  void SelectBead(int p);
  void Wave();
  void AddCircle(int nNano);
  void AddCylinder(int nNano);
  void AddPore(int nNano);
  void AddRigid();
  double HeightBoundary(double *Pos,int dir);
  //-----------ForceCreate.cpp----------------
  void CreateInitial();
  void  Create2d();
  void Create3d();
  void CreateLeaves();
  void CreateStalk();
  void CreatePore();
  void Create1d();
  void CreateRigid();
  void CreateMC();
  void CreateMD();
  void CreateRod();
  void CreateElectro();
  //---------ElPolyTens.cpp--------------------------------
  void AllocTens();
  void CalcTens();
  void CalcDens();
  double TensRefCart(double *Pos1,double *Pos2,double *PosP1,double *PosP2);
  double TensRefPol(double *Pos1,double *Pos2,double *PosP1,double *PosP2);
  typedef double(Forces::*TENS_REF)(double *Pos1,double *Pos2,double *PosP1,double *PosP2);
  TENS_REF Tens_Ref;
  double TensRef(double *Pos1,double *Pos2,double *PosP1,double *PosP2){
    return (*this.*Tens_Ref)(Pos1,Pos2,PosP1,PosP2);
  };
  void SumTens(int p1,int p2,double Forces,double *DistRel);
  void SumTens(int p1,int p2,double *Pre);
  void SumTens(double *Pos1,double *Pos2,double *Pre);
  void WriteTens(char *TFile,int Comp,double InvNFile);
  void WriteTens2d(FILE *FWrite,int Comp,double InvNFile);
  //---------ForcesLoop.cpp--------------------------------
  void ExplorePepSize();
  void ExplorePepSize2d();
  void ExploreDoubleMin();
  void CalcTotNrg(char *FName,int nFile);
  void CalcNrgPep(char *File2Open,int f);
  void RunDynamics();
  void RunWidom(char *File2Read,int f);
  void RunWidomChIn(char *File2Read,int f);
  void RunWidomChOut(char *File2Read,int f);
  void RosenIn(FILE *WidomIn);
  void RosenOut(FILE *WidomIn);
  void RunWidomBiasChOut(FILE *WidomOut);
  void ChooseSimMode();
  void Task();
  void CalcTens(char **argv,int *FilePos,int NFile);
  void AvForces(char **argv,int *FilePos,int NFile);
  void Trial();
  void MinimalMD();
  double MinimalNrg();
  void Sim1d(){VelVerlet1();ForceFieldLine();VelVerlet2();};
  void Sim2d(){VelVerlet1();Wave();VelVerlet2();};
  void Sim3d(){VelVerlet1();ForceFieldBulk();VelVerlet2();};
  void SimLeaves(){VelVerlet1();ForceFieldLeaves();AddRigid();VelVerlet2();};
  void SimRigid(){VelVerlet1();ForceFieldRigid();VelVerletRigid();VelVerlet2();};
  void MinimizeSol();
 //----------------CONSTANTS----------------------------
  double IncrDist;
  double Deltat;
  double NChemPotId;
  double ChemPotId;
  double ChemPotEx;
  double GaussVar;
  double CutOff;
  double Dx;
  double Viscosity;
  double Time;
  double NrgPBead;
  int Bead2Move;
  int Old2Move;
  int IntMax;
  int nEdge[3];
  int SysShape;
  int CalcMode;
  int ThermMode;
  int SimLimit;
  int SysAlloc;
  int IfInterp;
  int IfLeaves;
  int IfMove;
  int IfNano;
  int IfFillMatrix;
  int IfExit;
  int NInsertion;
  int NRemoval;
  int NFile[2];
  int NUpdate;
  int NWrite;
  int NSpline;
  KFORCES Kf;
  FORCES *Fm;
  TENS Tens;
  DomDec *Pc;
  PART *Pl;
#ifdef USE_MPI
  SingProc *Proc;
#endif
#ifdef __glut_h__

  int Interp();
  int Graphics(int argc,char **argv);
  void keyboard(unsigned char key,int x, int y);
  void DrawSoil();
  void DrawCarpet();
  void DrawParticles();
  void DrBondLine(int p);
  void DrawScene();
  void DrawNano();
  void DynamicsView();
  void Menu();
  GLuint *Cylinder;
  GLuint Particles;
  int NShow;
  int IfMovie;
  int Frame;
  int IfExt;
  int IfSpline;
  int IfRot;
  int IfSphere;
  int BeadType;
  int menu,submenu;
  int IfLine;
#endif
};
extern void DynamicsMotion();
#endif //ELPOLY_H