13   else if(VAR_IF_TYPE(
SysShape,SYS_LEAVES))
    15   else if(VAR_IF_TYPE(
SysShape,SYS_PORE))
    17   else if(VAR_IF_TYPE(
SysShape,SYS_ROD))
    27   else if(VAR_IF_TYPE(
SysShape,SYS_LEAVES)){
    30   else if(VAR_IF_TYPE(
SysShape,SYS_1D)){
    33   else if(VAR_IF_TYPE(
SysShape,SYS_3D)){
    36   else if(VAR_IF_TYPE(
SysShape,SYS_ROD)){
    39   else if(VAR_IF_TYPE(
SysShape,SYS_RIGID)){
    43   else if(VAR_IF_TYPE(
SysShape,SYS_MC)){
    48     else if(VAR_IF_TYPE(
CalcMode,CALC_NcVT)){
    51     else if(VAR_IF_TYPE(
CalcMode,CALC_mVT)){
    57     else if(VAR_IF_TYPE(
CalcMode,CALC_mcVT)){
    58       if(!VAR_IF_TYPE(
CalcMode,CALC_CONF_BIAS)){
    70   else if(VAR_IF_TYPE(
SysShape,SYS_ELECTRO)){
    74     else if(VAR_IF_TYPE(
CalcMode,CALC_mVT)){
    81   else if(VAR_IF_TYPE(
SysShape,SYS_MD)){
    84     for(
int p=0;p<
pNPart();p++){
    98   Shout(
"Explore pep size\n");
   103   double HeiStep = (HeiMax-HeiMin)/(
double)NGrid;
   106   double AngStep = (AngMax-AngMin)/(
double)NGrid;
   107   double SLapMin = 0.001;
   108   double SLapMax = 10.0;
   109   double SLapStep = 10.;
   110   for(
double SLap=SLapMin;SLap<SLapMax;SLap*=SLapStep){
   111     sprintf(FName,
"PepHeiElThickSLap%lf.dat",SLap);
   112     FILE *FWrite = fopen(FName,
"w");
   113     for(
double Hei=HeiMin;Hei<HeiMax;Hei+=HeiStep){
   114       for(
double Angle=AngMin;Angle<AngMax;Angle+=AngStep){
   115    fprintf(stderr,
"%lf %lf %lf\r",SLap,Hei,Angle);
   121    double Min = 
pEdge(2);
   126      if(
Pm[p].Pos[2] > Max && 
Pm[p].Typ == 0){
   133      if(
Pm[p].Pos[2] < Min && 
Pm[p].Typ == 0){
   138    fprintf(FWrite,
"%lf %lf %lf %lf %lf\n",Hei/.2,Angle,(Min-Max)/.2,SLap,.5*(xMin+xMax)/.2);
   146   Shout(
"Explore pep size 2d\n");
   152   double HeiStep = (HeiMax-HeiMin)/(
double)NGrid;
   155   double AngStep = (AngMax-AngMin)/(
double)NGrid;
   156   double InvNBin = 1./(double)NBin;
   158   for(
double Angle=AngMin;Angle<AngMax;Angle+=AngStep){
   159     for(
double Hei=HeiMin;Hei<HeiMax;Hei+=HeiStep){
   160       fprintf(stderr,
"%lf %lf\r",Hei,Angle);
   164       sprintf(FName,
"Sol2dHei%.1fAng%02d.dat",Hei,(
int)Angle);
   171   Shout(
"Explore pep size\n");
   181   FILE *FProf = fopen(
"DistThinMin.dat",
"w");
   182   for(
int i=0;i<NGrid;i++){
   183     fprintf(stderr,
"done %.2f %%\r",100.*i/(
double)NGrid);
   186     sprintf(FName,
"MinProf%03d.dat",i);
   187     FILE *FWrite = fopen(FName,
"w");
   189     fprintf(FWrite,String);
   191     fprintf(FWrite,String);
   192     double Min = 100000.;
   194     double zOld = 100000.;
   199       double z = 
pPos(p,2);
   202       fprintf(FWrite,
"%lf %lf\n",x,z);
   212     fprintf(FProf,
"%lf %lf %lf\n",
Nano[1].Pos[0]-
Nano[0].Pos[0],Dist,Min);
   215     sprintf(FName,
"Sol2dDist%.3f.dat",
Nano[1].Pos[0]-
Nano[0].Pos[0]);
   224   Shout(
"Widom on particles\n");
   228   sprintf(File2Write,
"WidomOut%05d.dat",f);
   229   ReOpen(File2Open,BF_PART);
   230   FILE *WidomOut = fopen(File2Write,
"w");
   231   SigErr(WidomOut==NULL,
"Cannot allocate %s\n",File2Write);
   232   for(
int p=0;p<
pNPart();p++){
   234     fprintf(WidomOut,
"%lf\n",NrgDiff);
   237   sprintf(File2Write,
"WidomIn%05d.dat",f);
   238   FILE *WidomIn = fopen(File2Write,
"w");
   239   for(
int p=0;p<NInt;p++){
   241     fprintf(WidomIn,
"%lf\n",NrgDiff);
   246   Shout(
"Widom insertion biased/Rosenbluth\n");
   249   for(
int p=0;p<NInt;p++){
   251     fprintf(WidomIn,
"%g\n",Weight);
   255   Shout(
"Widom deletion biased/Rosenbluth\n");
   260     fprintf(WidomIn,
"%g\n",Weight);
   264   Shout(
"Widom chain in\n");
   268   ReOpen(File2Open,BF_PART);
   269   sprintf(File2Write,
"WidomIn%05d.dat",f);
   270   FILE *WidomIn = fopen(File2Write,
"w");
   272   for(
int p=0;p<NInt;p++){
   275     fprintf(WidomIn ,
"%lf\n",NrgDiff[2]);
   280   Shout(
"Widom chain out\n");
   283   sprintf(File2Write,
"WidomOut%05d.dat",f);
   284   ReOpen(File2Open,BF_PART);
   285   FILE *WidomOut = fopen(File2Write,
"w");
   289     fprintf(WidomOut,
"%lf\n",NrgDiff[2]);
   294   Shout(
"Calculating total energy\n");
   295   ReOpen(File2Open,BF_PART);
   306       double CosAngle = 0.;
   307       for(
int d=0;d<3;d++){
   309    DistBA[d] -= floor(DistBA[d]/(
pEdge(d)) + .5)*
pEdge(d);
   311    DistCB[d] -= floor(DistCB[d]/(
pEdge(d)) + .5)*
pEdge(d);
   312    DistBA2 += SQR(DistBA[d]);
   313    DistCB2 += SQR(DistCB[d]);
   314    CosAngle += DistBA[d]*DistCB[d];
   316       DistCB2 = sqrt(DistCB2);
   318       if(p == c*
pNPCh()) 
continue;
   319       DistBA2 = sqrt(DistBA2);
   320       CosAngle /= (DistBA2*DistCB2);
   321       NrgBend += 
pkBen()*(1.-CosAngle);
   325   for(
int p=0;p<
pNPart();p++){
   329   printf(
"Nb %lf + Nano %lf = %lf Spr %lf + Ben %lf = %lf\n",OldNrgSys,NrgNano,OldNrgSys+NrgNano,NrgSpr,NrgBend,NrgSpr+NrgBend);
   330   fprintf(StatFile1,
"%d %lf %lf %lf %lf\n",f,OldNrgSys,NrgSpr,NrgBend,NrgNano);
   334   Shout(
"Calculating total energy\n");
   335   ReOpen(File2Open,BF_PART);
   340     if(!strncmp(
Block[b].Name,
"PEP",3)){
   344       for(
int p=pOff,link=0;p<MIN(pOff+
pNPCh(b),
pNPart());p++){
   350   printf(
"Nrg pep %lf\n",NrgPep);
   351   fprintf(StatFile1,
"%d %lf %lf\n",f,NrgPep,OldNrgSys);
   355   Shout(
"Calculating tension\n");
   358   for(
int f=0;f<
NFile;f++){
   359     fprintf(stderr,
"Elaborating file %s %.3f %%\r",argv[FilePos[f]],f/(
double)NFile*100.);
   360     ReOpen(argv[FilePos[f]],BF_PART);
   366      sprintf(FName,
"Tension2d%05dL%d.dat",f,c);
   368      sprintf(FName,
"Tension3d%05dL%d.dat",f,c);
   376   Shout(
"Calculating force field/average force\n");
   379   const int NBin = 120;
   380   double *Profile = (
double *)calloc(3*NBin,
sizeof(
double));
   381   for(
int f=0;f<
NFile;f++){
   382     fprintf(stderr,
"Elaborating file %s %.3f %%\r",argv[FilePos[f]],f/(
double)NFile*100.);
   383     ReOpen(argv[FilePos[f]],BF_PART);
   385     FILE *FForce = fopen(
"ForceProfile.dat",
"w");
   386     double Force[3] = {0.,0.,0.};
   388     for(
int p=0;p<
pNPart();p++){
   389       for(
int d=0;d<3;d++){
   390    cBin[d] = (int)(
Pm[p].Pos[d]*
pInvEdge(d)*NBin);
   391    if(cBin[d] < 0 || cBin[d] >= NBin) 
continue;
   392    Profile[cBin[d]*3+d] += 
Fm[p].
Dir[d];
   393    Force[d] += 
Fm[p].
Dir[d];
   394    fprintf(FForce,
"%d %g %g %g\n",p,Force[0],Force[1],Force[2]);
   397     printf(
"%g %g %g\n",Force[0],Force[1],Force[2]);
   400   FILE *FProf = fopen(
"ForceProfile.dat",
"w");
   401   for(
int b=0;b<NBin;b++){
   402     fprintf(FProf,
"%d %lf %lf %lf\n",b,Profile[b*3],Profile[b*3+1],Profile[b*3+2]);
   404   FILE *FField = fopen(
"ForceField.dat",
"w");
   405   double Delta = sqrt(
Kf.
CutOff2)/(double)NTab;
   406   for(
int b=0;b<NTab;b++){
   408     fprintf(FField,
"%lf ",x);
   409     for(
int t1=0;t1<
pNType();t1++){
   410       for(
int t2=0;t2<
pNType();t2++){
   411    fprintf(FField,
"%lf ",FTab[(b*
pNType()+t1)*
pNType()+t2]);
   414     fprintf(FField,
"\n");
   439   else if(VAR_IF_TYPE(
SysAlloc,ALL_MD)){
   442       for(
int p=0;p<
pNPart();p++){
   443    for(
int d=0;d<3;d++){
   444      v2 += SQR(
Pm[p].Vel[d]);
   456   Shout(
"run dynamics\n");
   457   double NLoopSec = 0.;
   465       CurrTime = time(NULL);
   466       NLoopSec += (s)/(
double)(CurrTime-InitTime);
   468       for(
int p=0;p<
pNPart();p++){
   469    for(
int d=0;d<3;d++){
   470      Temp += SQR(
Pm[p].Vel[d]);
   473       Temp = Temp/(double)(3*
pNPart());
   474       fprintf(stderr,
"NPart %d loop/ms %.3g acc/step %.3f in/out %.4f T %.4f accomplished %.3f %% Nrg %lf\n",
pNPart(),NLoopSec,(
NRemoval+
NInsertion)/(
double)s,
NInsertion/(
double)
NRemoval,Temp,s/(
double)(SimLimit)*100.,OldNrgSys);
   475       fprintf(StatFile1,
"%d %lf %lf\n",s,OldNrgSys,Temp);
   478       sprintf(FileName,
"Trajectory%09d.dat",
pStep());
   483   CurrTime = time(NULL);
   484   NLoopSec += (
SimLimit)/(
double)(CurrTime-InitTime);
   486   for(
int p=0;p<
pNPart();p++){
   487     for(
int d=0;d<3;d++){
   488       v2 += SQR(
Pm[p].Vel[d]);
   500   for(
int p1=0;p1<
pNPart();p1++){
   501     if(VAR_IF_TYPE(
SysAlloc,ALL_FORCES)){
   502       for(
int d=0;d<3;d++) 
Fm[p1].Dir[d] = 0.;
   506       if(p2 <= p1) 
continue;
   510       double InvDist = DistRel[3];
   511       if(VAR_IF_TYPE(
SysAlloc,ALL_FORCES)){
   512    for(
int d=0;d<3;d++){
   513      Fm[p1].
Dir[d] += Cons*DistRel[d]*InvDist;
   514      Fm[p2].
Dir[d] -= Cons*DistRel[d]*InvDist;
   522   double Sigma = sqrt(
pTemp());
   524   for(
int p=0;p<
pNPart();p++){
   525     for(
int d=0;d<3;d++){
   534   for(
int p=0;p<
pNPart();p++){
   535     for(
int d=0;d<3;d++){
   537       Temp += SQR(
Pm[p].Vel[d]);
   541   for(
int p=0;p<
pNPart();p++){
   543       for(
int d=0;d<3;d++){
   549   double Norm = 1./sqrt(Temp);
   557   double *Sol2d = (
double *)calloc(SQR(
NEdge),
sizeof(double));
   558   double *Count = (
double *)calloc(SQR(
NEdge),
sizeof(double));
   559   double NInvEdge = 1./(double)
NEdge;
   560   for(
int p=0;p<
pNPart();p++){
   561     if(
pType(p) != 0) 
continue;
   563     if(vx < 0 || vx >= 
NEdge) 
continue;
   565     if(vy < 0 || vy >= 
NEdge) 
continue;
   567     Count[vx*
NEdge+vx] = 1.;
   569   for(
int vx=0;vx<SQR(
NEdge);vx++){
   570     double Norm = Count[vx] > 0. ?  1./Count[vx] : 1.;
   575   for(
int n=0;n<NInt;n++){
 int TryMoveCh()
Trial desplacement of a chain. 
 
void AddPart(const int p, double *Pos)
Add a particle to the cell c. 
 
void RunDynamics()
Run a step further. 
 
void CalcTens()
Calculate the forces for the tension profile. 
 
void WidomBiasChIn(double *Weight)
Widom with Rosenbluth weight. 
 
double pkBen()
Bending coupling. 
 
void ExplorePepSize2d()
Find the minimun bilayer thickness for different peptide sizes. 
 
void CalcTotNrg(char *FName, int nFile)
Total energy of the system. 
 
double CalcTotNrgCh()
Calculate and sum up the energy of the chains. 
 
void VelVerlet1()
First step of the velocity Verlet. 
 
int NInsertion
Count accepted moves. 
 
int ForceFieldLine()
Helfrich Hamiltonian for a line. 
 
void WriteTens(char *TFile, int Comp, double InvNFile)
Write the pressure and density profile. 
 
void Create2d()
Create a plane of connected beads. 
 
BLOCK * Block
Information for every block. 
 
void Erase()
Erase the pairlist. 
 
double Vel[4]
xyzr Velocity of the particle 
 
int TryRemoveChBias()
Trial biased removal of a chain. 
 
NANO * Nano
Extra particle. 
 
double Hamaker
Strength of the interaction. 
 
void RunWidom(char *File2Read, int f)
Build the widom histograms. 
 
bool Write(char *OutFile)
Writes a "system-file" or a "x y z" file". 
 
void StringNano(char *NString, int n)
String for the rigid inclusion in the header file. 
 
void CheckPairList()
Check the pair list. 
 
void SetTime(double Val)
Set Time. 
 
double pInvEdge(int d)
Inverted xyzr edges of the simulation box. 
 
double Height
Height of the cylinder. 
 
int pType(int p)
Return the type. 
 
void Dynamics()
Sum up all the forces and update the positions. 
 
int TryInsertCh()
Trial insertion of a chain. 
 
double Pos[3]
xyz Position of the particle 
 
void RunWidomChIn(char *File2Read, int f)
Build the widom histograms. 
 
double SumForcesMD()
Iterate all over the particles and calculate the forces. 
 
void WidomInsert(double *NrgDiff)
Widom insertion. 
 
int NWrite
How many timesteps before write the snapshot. 
 
double pkSpr()
Spring coupling. 
 
double pSprRest()
Rest distance of the harmonic potential. 
 
void CalcNrgPep(char *File2Open, int f)
Exchange energy of the protein. 
 
void ForceFieldRod()
Bending potential on a rod. 
 
void Dist2Curr(double *DistRel)
Retrun the squared current interparticle distance. 
 
void SolveLinksIterative()
Solve a system of four oder differential equation of particles connected by links. 
 
int TryInsertChBias()
Trial biased insertion of a chain. 
 
void WidomInsertCh(double *NrgDiff)
Widom insertion. 
 
int NFile[2]
First and last file of the list. 
 
int SysAlloc
Which arrays are allocated. 
 
void WidomRemoveCh(double *NrgDiff, int c)
Widom removal. 
 
int IfFillMatrix
If the matrix has to be changed. 
 
void CalcForcesDensFunc()
Calculate the forces for the density functional. 
 
int NRemoval
Count accepted moves. 
 
int SimLimit
Maximum number of time steps. 
 
void ExploreDoubleMin()
Find the minimum for different interpeptide distances. 
 
int CalcMode
Calculation mode. 
 
void ApplyTherm()
Pointer to the energy function. 
 
void Shout(const char *s,...)
Internal message. 
 
double pNanoPos(int n, int d)
Return back folded nano position. 
 
double pEdge(int d)
xyzr edges of the simulation box 
 
void VelVerletRigid()
Velocity Verlet for a rigid body, first step. 
 
int pNPCh()
Number of particle per chain. 
 
void ForceFieldRigid()
Interaction between rigid bodies. 
 
int pNBlock()
Number of blocks. 
 
void RosenOut(FILE *WidomIn)
Rosenbluth histograms for deletion. 
 
double NanoNrg(int p)
Exchange energy with the nano. 
 
int TryMove()
Trial movement. 
 
void ReOpen(char *FName, int Bf)
Open a new file. 
 
int NEdge
Number of particles per edge. 
 
double Casuale()
Random uniform number. 
 
bool WriteTxvl(char *OutFile)
Writes a "system-file" or a "x y z" file". 
 
KFORCES Kf
Prefactor of the forces. 
 
int NUpdate
How many timesteps before redrawing. 
 
void NextCurr()
Increase the iterator to the next couple. 
 
double MinimalNrg()
Minmal nrg. 
 
int ForceFieldLeaves()
Helfrich Hamiltonian with an elastic coupling. 
 
int TryRemove()
Trial removal. 
 
void ExplorePepSize()
Find the minimun bilayer thickness for different peptide sizes. 
 
int SysShape
Shape of system. 
 
void RosenIn(FILE *WidomIn)
Rosenbluth weights for insertion. 
 
int pNNano()
Number of nanoparticles. 
 
void RunWidomChOut(char *File2Read, int f)
Build the widom histograms. 
 
void MinimalMD()
Minmal md. 
 
int TryRemoveCh()
Trial removal of a chain. 
 
int nEdge[3]
Number of particle per edge. 
 
void IncrStep()
Increment Step. 
 
int ForceFieldBulk()
Armonic potential on a lattice. 
 
void MinimizeSol()
Iterative process to approach to the solution. 
 
double Potential(double Dist, int t1, int t2, double *Pot)
Pointer to a potential. 
 
void SolveLeaves()
Solve a system of four oder differential equation of particles in a line. 
 
int CalcMode
Calculation mode. 
 
void ClearDens()
Set the local densities to zero. 
 
double Gaussiano(double Media, double Scarto)
Gaussian random number. 
 
int SysFormat
Contains the definition of the file format. 
 
void Wave()
Sinusoidal surface wave. 
 
double CalcNrgCh(int c, double *Pot)
Pointer to the chain energy function. 
 
void AllocTens()
Alloc the pressure profile. 
 
void VelVerlet2()
Second step of the velocity Verlet. 
 
Matematica * Mat
Implementation of all usefull algorythms. 
 
int pNType()
of types of the particle
 
void Task()
Perform a operation every time step. 
 
double pPos(int p, int d)
Return back folded position. 
 
double pTemp()
Temperature. 
 
int pNChain()
Number of chain. 
 
int IfCurr()
Tell when the curr loop is over. 
 
void WidomRemove(double *NrgDiff, int p)
Widom removal. 
 
void SolveRod()
Solve a system of four oder differential equation of particles connected by links. 
 
double CutOff2
CutOff of the lennard jones. 
 
int p2Curr
Current particle. 
 
double SumDens(int pInit, int pEnd)
Sum the local density for the particles between pInit and pEnd and multiply the factors by the virial...
 
void ChooseSimMode()
Choose the simulation method. 
 
void CalcDens()
Calculate the densities. 
 
PART * Pm
Particle information of all particle. 
 
int AddDens(int pInit, int pEnd)
Add the densities connected with the particles between pInit and pEnd. 
 
void SetCurr(int p)
Gather information of the neighbouring cells. 
 
TENS Tens
Structure for the pressure calculation. 
 
double SLap
Prefactor of the square laplacian. 
 
void AvForces(char **argv, int *FilePos, int NFile)
Average of the forces. 
 
void SetEdge(double Val, int d)
Set Edge. 
 
int TryInsert()
Trial insertion. 
 
FORCES * Fm
Array containing the forces for each particle. 
 
int pStep()
Number of steps. 
 
int pNPart()
Number of particle. 
 
void DefForceParam()
Define the parameters for calculating the force. 
 
void WidomBiasChOut(double *Weight, int c)
Widom with Rosenbluth weight. 
 
void Solve()
Solve a system of four oder differential equation.