X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=moldyn.c;h=48129029a1c6e47a6fc0a59893c9fef04caa8999;hb=c2697c8b6f808b7edef2f9e2ceded0aa65374168;hp=82ae543c073457a47372253f6c238f76d1bfe6b5;hpb=0d1dfb1e5027d215fced8ca306dd658f692a2a44;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 82ae543..4812902 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -30,101 +30,6 @@
 #include "potentials/tersoff.h"
 #endif
 
-
-/*
- * global variables, pse and atom colors (only needed here)
- */ 
-
-static char *pse_name[]={
-	"*",
-	"H",
-	"He",
-	"Li",
-	"Be",
-	"B",
-	"C",
-	"N",
-	"O",
-	"F",
-	"Ne",
-	"Na",
-	"Mg",
-	"Al",
-	"Si",
-	"P",
-	"S",
-	"Cl",
-	"Ar",
-};
-
-static char *pse_col[]={
-	"*",
-	"White",
-	"He",
-	"Li",
-	"Be",
-	"B",
-	"Gray",
-	"N",
-	"Blue",
-	"F",
-	"Ne",
-	"Na",
-	"Mg",
-	"Al",
-	"Yellow",
-	"P",
-	"S",
-	"Cl",
-	"Ar",
-};
-
-/*
-static double pse_mass[]={
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	M_C,
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	M_SI,
-	0,
-	0,
-	0,
-	0,
-};
-
-static double pse_lc[]={
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	LC_C,
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	0,
-	LC_SI,
-	0,
-	0,
-	0,
-	0,
-};
-*/
-
 /*
  * the moldyn functions
  */
@@ -208,7 +113,7 @@ int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) {
 	moldyn->pt_scale|=ptype;
 	moldyn->p_tc=ptc;
 
-	printf("[moldyn] p/t scaling:\n");
+	printf("[moldyn] p scaling:\n");
 
 	printf("  p: %s",ptype?"yes":"no ");
 	if(ptype)
@@ -224,7 +129,7 @@ int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) {
 	moldyn->pt_scale|=ttype;
 	moldyn->t_tc=ttc;
 
-	printf("[moldyn] p/t scaling:\n");
+	printf("[moldyn] t scaling:\n");
 
 	printf("  t: %s",ttype?"yes":"no ");
 	if(ttype)
@@ -315,7 +220,7 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 			moldyn->func3b_k1=tersoff_mult_3bp_k1;
 			moldyn->func3b_j2=tersoff_mult_3bp_j2;
 			moldyn->func3b_k2=tersoff_mult_3bp_k2;
-			// missing: check 2b bond func
+			moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
 			break;
 		case MOLDYN_POTENTIAL_AM:
 			moldyn->func3b_j1=albe_mult_3bp_j1;
@@ -580,6 +485,10 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	count=moldyn->count;
 
 	/* how many atoms do we expect */
+	if(type==NONE) {
+		new*=1;
+		printf("[moldyn] WARNING: create 'none' lattice called");
+	}
 	if(type==CUBIC) new*=1;
 	if(type==FCC) new*=4;
 	if(type==DIAMOND) new*=8;
@@ -963,7 +872,7 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 		scale*=2.0;
 	else
 		if(moldyn->pt_scale&T_SCALE_BERENDSEN)
-			scale=1.0+(scale-1.0)/moldyn->t_tc;
+			scale=1.0+(scale-1.0)*moldyn->timestep/moldyn->t_tc;
 	scale=sqrt(scale);
 
 	/* velocity scaling */
@@ -1273,8 +1182,8 @@ int scale_volume(t_moldyn *moldyn) {
 
 	/* scaling factor */
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
-		scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc;
-		scale=pow(scale,ONE_THIRD);
+		scale=(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->timestep;
+		scale=pow(1.0-scale,ONE_THIRD);
 	}
 	else {
 		scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
@@ -1669,13 +1578,16 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
 	if(moldyn->cutoff>0.5*moldyn->dim.z)
 		printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
-	ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
-	if(ds>0.05*moldyn->nnd)
+	if(moldyn->count) {
+		ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
+		if(ds>0.05*moldyn->nnd)
 		printf("[moldyn] WARNING: forces too high / tau too small!\n");
+	}
 
 	/* zero absolute time */
-	moldyn->time=0.0;
-	moldyn->total_steps=0;
+	// should have right values!
+	//moldyn->time=0.0;
+	//moldyn->total_steps=0;
 
 	/* debugging, ignore */
 	moldyn->debug=0;
@@ -1792,6 +1704,7 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
        moldyn->p/BAR,moldyn->p_avg/BAR,
+       //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));
 
@@ -2212,7 +2125,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 		/* check forces regarding the given timestep */
 		if(v3_norm(&(itom[i].f))>\
-		   0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
+                   0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
 			printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
 			       i);
 	}
@@ -2681,8 +2594,6 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 	if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
 		return 0;
 
-	d=sqrt(d);
-
 	/* now count this bonding ... */
 	ba=data;