changed energy fluctuation calc and heat cap calc

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index e554a2f..cd2a803 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -404,11 +404,14 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
        if(moldyn->rfd) {
                dprintf(moldyn->rfd,report_end);
                close(moldyn->rfd);
-               snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
+               snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
+                        moldyn->vlsdir);
                system(sc);
-               snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
+               snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
+                        moldyn->vlsdir);
                system(sc);
-               snprintf(sc,255,"cd %s && dvipdf report",moldyn->vlsdir);
+               snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1",
+                        moldyn->vlsdir);
                system(sc);
        }
        if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
@@ -501,6 +504,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
                check_per_bound(moldyn,&(atom[ret].r));
        }
 
+       /* update total system mass */
+       total_mass_calc(moldyn);
+
        return ret;
 }
 
@@ -642,6 +648,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        atom[count].tag=count;
        atom[count].attr=attr;
 
+       /* update total system mass */
+       total_mass_calc(moldyn);
+
        return 0;
 }
 
@@ -702,6 +711,18 @@ int thermal_init(t_moldyn *moldyn,u8 equi_init) {
        return 0;
 }
 
+double total_mass_calc(t_moldyn *moldyn) {
+
+       int i;
+
+       moldyn->mass=0.0;
+
+       for(i=0;i<moldyn->count;i++)
+               moldyn->mass+=moldyn->atom[i].mass;
+
+       return moldyn->mass;
+}
+
 double temperature_calc(t_moldyn *moldyn) {
 
        /* assume up to date kinetic energy, which is 3/2 N k_B T */
@@ -822,7 +843,56 @@ double pressure_calc(t_moldyn *moldyn) {
        moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps;
 
        return moldyn->p;
-}      
+}
+
+int energy_fluctuation_calc(t_moldyn *moldyn) {
+
+       /* assume up to date energies */
+
+       /* kinetic energy */
+       moldyn->k_sum+=moldyn->ekin;
+       moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
+       moldyn->k_mean=moldyn->k_sum/moldyn->total_steps;
+       moldyn->k2_mean=moldyn->k2_sum/moldyn->total_steps;
+       moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean);
+
+       /* potential energy */
+       moldyn->v_sum+=moldyn->energy;
+       moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
+       moldyn->v_mean=moldyn->v_sum/moldyn->total_steps;
+       moldyn->v2_mean=moldyn->v2_sum/moldyn->total_steps;
+       moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean);
+
+       return 0;
+}
+
+int get_heat_capacity(t_moldyn *moldyn) {
+
+       double temp2,ighc;
+
+       /* (temperature average)^2 */
+       temp2=moldyn->mean_t*moldyn->mean_t;
+       printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
+              moldyn->mean_t);
+
+       /* ideal gas contribution */
+       ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
+       printf("  ideal gas contribution: %f\n",
+              ighc/moldyn->mass*KILOGRAM/JOULE);
+
+       /* specific heat for nvt ensemble */
+       moldyn->c_v_nvt=moldyn->dv2_mean/(K_BOLTZMANN*temp2)+ighc;
+       moldyn->c_v_nvt/=moldyn->mass;
+
+       /* specific heat for nve ensemble */
+       moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_mean/(ighc*K_BOLTZMANN*temp2)));
+       moldyn->c_v_nve/=moldyn->mass;
+
+       printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
+       printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
+
+       return 0;
+}
 
 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
@@ -1304,6 +1374,7 @@ return 0;
                e_kin_calc(moldyn);
                temperature_calc(moldyn);
                pressure_calc(moldyn);
+               energy_fluctuation_calc(moldyn);
                //tp=thermodynamic_pressure_calc(moldyn);
 //printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR);
 
@@ -1376,6 +1447,8 @@ return 0;
                               moldyn->mean_gp/BAR,
                               moldyn->volume);
                        fflush(stdout);
+printf("\n");
+get_heat_capacity(moldyn);
                }
 
                /* increase absolute time */
@@ -1385,8 +1458,9 @@ return 0;
        }
 
                /* check for hooks */
-               if(sched->hook)
-                       sched->hook(moldyn,sched->hook_params);
+               if(sched->count+1<sched->total_sched)
+                       if(sched->hook)
+                               sched->hook(moldyn,sched->hook_params);
 
                /* get a new info line */
                printf("\n");
@@ -1805,3 +1879,26 @@ int moldyn_bc_check(t_moldyn *moldyn) {
 
        return 0;
 }
+
+/*
+ * post processing functions
+ */
+
+int get_line(int fd,char *line,int max) {
+
+       int count,ret;
+
+       count=0;
+
+       while(1) {
+               if(count==max) return count;
+               ret=read(fd,line+count,1);
+               if(ret<=0) return ret;
+               if(line[count]=='\n') {
+                       line[count]='\0';
+                       return count+1;
+               }
+               count+=1;
+       }
+}
+