trisurf-ng.git - Gitblit

Samo Penic

2019-03-08 2afc2f4f1dd89518995f1b5a539aea932aecab65

commit \| author \| age
7f6076	1	/* vim: set ts=4 sts=4 sw=4 noet : */
aec47d	2	#include<stdlib.h>
SP	3	#include<math.h>
	4	#include "general.h"
	5	#include "vertex.h"
	6	#include "bond.h"
	7	#include "triangle.h"
	8	#include "vesicle.h"
	9	#include "energy.h"
	10	#include "timestep.h"
	11	#include "cell.h"
9166cb	12	#include "io.h"
aec47d	13	#include<stdio.h>
SP	14	#include "vertexmove.h"
1ad6d1	15	#include <string.h>
43c042	16	#include "constvol.h"
51b4f0	17	#include "plugins.h"
aec47d	18
fedf2b	19	ts_bool single_verticle_timestep(ts_vesicle vesicle,ts_vertex vtx,ts_double *rn){
aec47d	20	ts_uint i;
SP	21	ts_bool retval;
	22	ts_uint cellidx;
77a2c7	23	ts_double delta_energy, oenergy,dvol=0.0, darea=0.0, dstretchenergy=0.0;
ed31fe	24	ts_double costheta,sintheta,phi,r;
1ad6d1	25	//This will hold all the information of vtx and its neighbours
2afc2f	26	ts_vertex backupvtx[20]; // constvol_vtx_moved=NULL, constvol_vtx_backup=NULL;
b5cd8c	27	memcpy((void )&backupvtx[0],(void )vtx,sizeof(ts_vertex));
672ae4	28
b5cd8c	29	//random move in a sphere with radius stepsize:
SP	30	r=vesicle->stepsize*rn[0];
	31	phi=rn[1]2M_PI;
	32	costheta=2*rn[2]-1;
	33	sintheta=sqrt(1-pow(costheta,2));
	34	vtx->x=vtx->x+rsinthetacos(phi);
	35	vtx->y=vtx->y+rsinthetasin(phi);
	36	vtx->z=vtx->z+r*costheta;
672ae4	37
fe24d2	38
51b4f0	39	/* Entry point for plugin vm_hard_constraint() function */
77a2c7	40	vesicle->plist->pointer=vesicle->plist->chain->vm_hard_constraint;
SP	41	while(vesicle->plist->pointer!=NULL){
	42	retval = vesicle->plist->pointer->plugin->function->vm_hard_constraint(vesicle,vtx, &backupvtx[0]);
51b4f0	43	if(retval==TS_FAIL){
SP	44	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	45	return TS_FAIL;
	46	}
77a2c7	47	vesicle->plist->pointer=vesicle->plist->pointer->next;
51b4f0	48	}
36dba8	49	/* End of vm_hard_constraint() */
51b4f0	50
36dba8	51	/* Backuping the neighbours */
1ad6d1	52	for(i=0;i<vtx->neigh_no;i++){
dcd350	53	memcpy((void )&backupvtx[i+1],(void )vtx->neigh[i],sizeof(ts_vertex));
1ad6d1	54	}
aec47d	55
36dba8	56	/* Entry point for plugin vm_energy_before_prepare() */
SP	57
	58	vesicle->plist->pointer=vesicle->plist->chain->vm_energy_before_prepare;
	59	while(vesicle->plist->pointer!=NULL){
	60	vesicle->plist->pointer->plugin->function->vm_energy_before_prepare(vesicle, vtx);
	61	vesicle->plist->pointer=vesicle->plist->pointer->next;
c0ae90	62	}
SP	63
	64	if(vesicle->tape->constareaswitch==2){
	65	for(i=0;i<vtx->tristar_no;i++) darea-=vtx->tristar[i]->area;
	66
	67	}
7ec6fb	68	//stretching energy 1 of 3
SP	69	if(vesicle->tape->stretchswitch==1){
699ac4	70	for(i=0;i<vtx->tristar_no;i++) dstretchenergy-=vtx->tristar[i]->energy;
7ec6fb	71	}
aec47d	72	delta_energy=0;
0dd5ba	73
SP	74
fe5069	75	// vesicle_volume(vesicle);
SP	76	// fprintf(stderr,"Volume in the beginning=%1.16e\n", vesicle->volume);
43c042	77
aec47d	78	//update the normals of triangles that share bead i.
8f6a69	79	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
a63f17	80	oenergy=vtx->energy;
aec47d	81	energy_vertex(vtx);
a63f17	82	delta_energy=vtx->xk*(vtx->energy - oenergy);
aec47d	83	//the same is done for neighbouring vertices
8f6a69	84	for(i=0;i<vtx->neigh_no;i++){
SP	85	oenergy=vtx->neigh[i]->energy;
	86	energy_vertex(vtx->neigh[i]);
	87	delta_energy+=vtx->neigh[i]->xk*(vtx->neigh[i]->energy-oenergy);
aec47d	88	}
414b8a	89
77a2c7	90
SP	91	/* Entry point for plugin vm_energy_after_execute() */
	92
	93	vesicle->plist->pointer=vesicle->plist->chain->vm_energy_after_execute;
	94	while(vesicle->plist->pointer!=NULL){
	95	delta_energy+=vesicle->plist->pointer->plugin->function->vm_energy_after_execute(vesicle, vtx);
	96	vesicle->plist->pointer=vesicle->plist->pointer->next;
	97	}
	98
	99
43c042	100
c0ae90	101	if(vesicle->tape->constareaswitch==2){
SP	102	/* check whether the darea is gt epsarea */
	103	for(i=0;i<vtx->tristar_no;i++) darea+=vtx->tristar[i]->area;
	104	if(fabs(vesicle->area+darea-A0)>epsarea){
	105	//restore old state.
	106	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	107	for(i=0;i<vtx->neigh_no;i++){
	108	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
	109	}
	110	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
	111	//fprintf(stderr,"fajlam!\n");
	112	return TS_FAIL;
	113	}
	114
	115
	116	}
1121fa	117
250de4	118	/* Vertices with spontaneous curvature may have spontaneous force perpendicular to the surface of the vesicle. additional delta energy is calculated in this function */
SP	119	delta_energy+=direct_force_energy(vesicle,vtx,backupvtx);
7ec6fb	120
SP	121	//stretching energy 2 of 3
	122	if(vesicle->tape->stretchswitch==1){
	123	for(i=0;i<vtx->tristar_no;i++){
	124	stretchenergy(vesicle, vtx->tristar[i]);
699ac4	125	dstretchenergy+=vtx->tristar[i]->energy;
7ec6fb	126	}
SP	127	}
	128
	129	delta_energy+=dstretchenergy;
	130
304510	131	/* No poly-bond energy for now!
fedf2b	132	if(vtx->grafted_poly!=NULL){
M	133	delta_energy+=
	134	(pow(sqrt(vtx_distance_sq(vtx, vtx->grafted_poly->vlist->vtx[0])-1),2)-
	135	pow(sqrt(vtx_distance_sq(&backupvtx[0], vtx->grafted_poly->vlist->vtx[0])-1),2)) *vtx->grafted_poly->k;
	136	}
304510	137	*/
0dd5ba	138
SP	139	// plane confinement energy due to compressing force
	140	if(vesicle->tape->plane_confinement_switch){
	141	if(vesicle->confinement_plane.force_switch){
	142	//substract old energy
	143	if(abs(vesicle->tape->plane_d/2.0-vesicle->confinement_plane.z_max)>1e-10) {
	144	delta_energy-=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_max-backupvtx[0].z,2);
	145	delta_energy+=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_max-vtx->z,2);
	146	}
	147	if(abs(-vesicle->tape->plane_d/2.0-vesicle->confinement_plane.z_min)>1e-10) {
	148	delta_energy-=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_min-backupvtx[0].z,2);
	149	delta_energy+=vesicle->tape->plane_F / pow(vesicle->confinement_plane.z_min-vtx->z,2);
	150	}
	151	}
	152	}
	153
77a2c7	154
SP	155	/* Entry point for plugin vm_before_montecarlo_constraint() function */
	156	vesicle->plist->pointer=vesicle->plist->chain->vm_before_montecarlo_constraint;
	157	while(vesicle->plist->pointer!=NULL){
	158	retval = vesicle->plist->pointer->plugin->function->vm_before_montecarlo_constraint(vesicle,vtx, &backupvtx[0]);
	159	if(retval==TS_FAIL){
	160	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
	161	for(i=0;i<vtx->neigh_no;i++){
	162	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
	163	}
	164	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
	165	return TS_FAIL;
	166	}
	167	vesicle->plist->pointer=vesicle->plist->pointer->next;
	168	}
	169	/* End of vm_before_montecarlo_constraint() */
	170
	171
	172
	173
	174
314f2d	175	// fprintf(stderr, "DE=%f\n",delta_energy);
aec47d	176	//MONTE CARLOOOOOOOO
SP	177	if(delta_energy>=0){
	178	#ifdef TS_DOUBLE_DOUBLE
3de289	179	if(exp(-delta_energy)< drand48())
aec47d	180	#endif
SP	181	#ifdef TS_DOUBLE_FLOAT
	182	if(expf(-delta_energy)< (ts_float)drand48())
	183	#endif
	184	#ifdef TS_DOUBLE_LONGDOUBLE
	185	if(expl(-delta_energy)< (ts_ldouble)drand48())
	186	#endif
	187	{
	188	//not accepted, reverting changes
fbcbdf	189	// fprintf(stderr,"MC failed\n");
dcd350	190	vtx=memcpy((void )vtx,(void )&backupvtx[0],sizeof(ts_vertex));
1ad6d1	191	for(i=0;i<vtx->neigh_no;i++){
a63f17	192	vtx->neigh[i]=memcpy((void )vtx->neigh[i],(void )&backupvtx[i+1],sizeof(ts_vertex));
1ad6d1	193	}
SP	194
aec47d	195	//update the normals of triangles that share bead i.
dcd350	196	for(i=0;i<vtx->tristar_no;i++) triangle_normal_vector(vtx->tristar[i]);
7ec6fb	197	//stretching energy 3 of 3
SP	198	if(vesicle->tape->stretchswitch==1){
	199	for(i=0;i<vtx->tristar_no;i++){
	200	stretchenergy(vesicle,vtx->tristar[i]);
	201	}
	202	}
1ad6d1	203
fe5069	204	// fprintf(stderr, "before vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
2afc2f	205	// if(vesicle->tape->constvolswitch == 1){
SP	206	// constvolumerestore(constvol_vtx_moved,constvol_vtx_backup);
	207	// }
fe5069	208	// fprintf(stderr, "after vtx(x,y,z)=%e,%e,%e\n",constvol_vtx_moved->x, constvol_vtx_moved->y, constvol_vtx_moved->z);
dd5aca	209	// vesicle_volume(vesicle);
SP	210	// fprintf(stderr,"Volume after fail=%1.16e\n", vesicle->volume);
aec47d	211	return TS_FAIL;
SP	212	}
	213	}
2b14da	214	//accepted
fbcbdf	215	// fprintf(stderr,"MC accepted\n");
a63f17	216	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
51b4f0	217	cellidx=vertex_self_avoidance(vesicle, vtx);
a63f17	218	if(vtx->cell!=vesicle->clist->cell[cellidx]){
SP	219	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	220	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	221	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx[0].cell,vtx);
	222
	223	}
2b14da	224
1121fa	225	if(vesicle->tape->constvolswitch == 2){
SP	226	vesicle->volume+=dvol;
2afc2f	227	}
SP	228
	229
c0ae90	230
SP	231	if(vesicle->tape->constareaswitch==2){
	232	vesicle->area+=darea;
	233	}
2afc2f	234
SP	235
	236	/* Entry point for plugin vm_before_montecarlo_constraint() function */
	237	vesicle->plist->pointer=vesicle->plist->chain->vm_new_state_accepted;
	238	while(vesicle->plist->pointer!=NULL){
	239	vesicle->plist->pointer->plugin->function->vm_new_state_accepted(vesicle,vtx, &backupvtx[0]);
	240	vesicle->plist->pointer=vesicle->plist->pointer->next;
	241	}
	242	/* End of vm_before_montecarlo_constraint() */
	243
	244
	245
	246
	247
aec47d	248	return TS_SUCCESS;
SP	249	}
	250
fedf2b	251
M	252	ts_bool single_poly_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	253	ts_uint i;
	254	ts_bool retval;
	255	ts_uint cellidx;
304510	256	// ts_double delta_energy;
fedf2b	257	ts_double costheta,sintheta,phi,r;
304510	258	ts_double dist;
fedf2b	259	//This will hold all the information of vtx and its neighbours
M	260	ts_vertex backupvtx;
304510	261	// ts_bond backupbond[2];
fedf2b	262	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	263
	264	//random move in a sphere with radius stepsize:
	265	r=vesicle->stepsize*rn[0];
	266	phi=rn[1]2M_PI;
	267	costheta=2*rn[2]-1;
	268	sintheta=sqrt(1-pow(costheta,2));
	269	vtx->x=vtx->x+rsinthetacos(phi);
	270	vtx->y=vtx->y+rsinthetasin(phi);
	271	vtx->z=vtx->z+r*costheta;
	272
	273
	274	//distance with neighbours check
304510	275	for(i=0;i<vtx->neigh_no;i++){
M	276	dist=vtx_distance_sq(vtx,vtx->neigh[i]);
	277	if(dist<1.0 \|\| dist>vesicle->dmax) {
	278	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	279	return TS_FAIL;
	280	}
	281	}
	282
	283	// Distance with grafted vesicle-vertex check:
	284	if(vtx==poly->vlist->vtx[0]){
	285	dist=vtx_distance_sq(vtx,poly->grafted_vtx);
	286	if(dist<1.0 \|\| dist>vesicle->dmax) {
	287	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	288	return TS_FAIL;
	289	}
	290	}
	291
fedf2b	292
M	293	//self avoidance check with distant vertices
	294	cellidx=vertex_self_avoidance(vesicle, vtx);
	295	//check occupation number
	296	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	297
	298	if(retval==TS_FAIL){
	299	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	300	return TS_FAIL;
	301	}
	302
	303
	304	//if all the tests are successful, then energy for vtx and neighbours is calculated
304510	305	/* Energy ignored for now!
fedf2b	306	delta_energy=0;
M	307	for(i=0;i<vtx->bond_no;i++){
	308	memcpy((void )&backupbond[i],(void )vtx->bond[i],sizeof(ts_bond));
	309
	310	vtx->bond[i]->bond_length=sqrt(vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2));
	311	bond_energy(vtx->bond[i],poly);
	312	delta_energy+= vtx->bond[i]->energy - backupbond[i].energy;
	313	}
	314
	315	if(vtx==poly->vlist->vtx[0]){
	316	delta_energy+=
	317	(pow(sqrt(vtx_distance_sq(vtx, poly->grafted_vtx)-1),2)-
	318	pow(sqrt(vtx_distance_sq(&backupvtx, poly->grafted_vtx)-1),2)) *poly->k;
	319
	320	}
	321
	322
	323	if(delta_energy>=0){
	324	#ifdef TS_DOUBLE_DOUBLE
	325	if(exp(-delta_energy)< drand48() )
	326	#endif
	327	#ifdef TS_DOUBLE_FLOAT
	328	if(expf(-delta_energy)< (ts_float)drand48())
	329	#endif
	330	#ifdef TS_DOUBLE_LONGDOUBLE
	331	if(expl(-delta_energy)< (ts_ldouble)drand48())
	332	#endif
	333	{
	334	//not accepted, reverting changes
	335	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	336	for(i=0;i<vtx->bond_no;i++){
	337	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
	338	}
	339
	340	return TS_FAIL;
	341	}
	342	}
304510	343	*/
fedf2b	344
M	345	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
	346	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	347	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	348	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	349	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	350	}
	351	// if(oldcellidx);
	352	//END MONTE CARLOOOOOOO
	353	return TS_SUCCESS;
	354	}
58230a	355
M	356
	357
	358
	359	ts_bool single_filament_vertex_move(ts_vesicle vesicle,ts_poly poly,ts_vertex vtx,ts_double rn){
	360	ts_uint i;
	361	ts_bool retval;
	362	ts_uint cellidx;
b30f45	363	ts_double delta_energy;
58230a	364	ts_double costheta,sintheta,phi,r;
M	365	ts_double dist[2];
	366	//This will hold all the information of vtx and its neighbours
b30f45	367	ts_vertex backupvtx,backupneigh[2];
58230a	368	ts_bond backupbond[2];
b30f45	369
M	370	//backup vertex:
58230a	371	memcpy((void )&backupvtx,(void )vtx,sizeof(ts_vertex));
M	372
	373	//random move in a sphere with radius stepsize:
	374	r=vesicle->stepsize*rn[0];
	375	phi=rn[1]2M_PI;
	376	costheta=2*rn[2]-1;
	377	sintheta=sqrt(1-pow(costheta,2));
	378	vtx->x=vtx->x+rsinthetacos(phi);
	379	vtx->y=vtx->y+rsinthetasin(phi);
	380	vtx->z=vtx->z+r*costheta;
	381
	382
	383	//distance with neighbours check
	384	for(i=0;i<vtx->bond_no;i++){
	385	dist[i]=vtx_distance_sq(vtx->bond[i]->vtx1,vtx->bond[i]->vtx2);
	386	if(dist[i]<1.0 \|\| dist[i]>vesicle->dmax) {
	387	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	388	return TS_FAIL;
	389	}
	390	}
	391
fe24d2	392	// TODO: Maybe faster if checks only nucleus-neighboring cells
M	393	// Nucleus penetration check:
	394	if (vtx->xvtx->x + vtx->yvtx->y + vtx->z*vtx->z < vesicle->R_nucleus){
	395	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	396	return TS_FAIL;
	397	}
	398
58230a	399
M	400	//self avoidance check with distant vertices
	401	cellidx=vertex_self_avoidance(vesicle, vtx);
	402	//check occupation number
	403	retval=cell_occupation_number_and_internal_proximity(vesicle->clist,cellidx,vtx);
	404	if(retval==TS_FAIL){
	405	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
	406	return TS_FAIL;
	407	}
	408
	409	//backup bonds
	410	for(i=0;i<vtx->bond_no;i++){
	411	memcpy(&backupbond[i],vtx->bond[i], sizeof(ts_bond));
	412	vtx->bond[i]->bond_length=sqrt(dist[i]);
	413	bond_vector(vtx->bond[i]);
b30f45	414	}
M	415
	416	//backup neighboring vertices:
	417	for(i=0;i<vtx->neigh_no;i++){
	418	memcpy(&backupneigh[i],vtx->neigh[i], sizeof(ts_vertex));
58230a	419	}
M	420
	421	//if all the tests are successful, then energy for vtx and neighbours is calculated
b30f45	422	delta_energy=0;
M	423
	424	if(vtx->bond_no == 2){
	425	vtx->energy = -(vtx->bond[0]->xvtx->bond[1]->x + vtx->bond[0]->yvtx->bond[1]->y + vtx->bond[0]->z*vtx->bond[1]->z)/vtx->bond[0]->bond_length/vtx->bond[1]->bond_length;
	426	delta_energy += vtx->energy - backupvtx.energy;
58230a	427	}
M	428
b30f45	429	for(i=0;i<vtx->neigh_no;i++){
M	430	if(vtx->neigh[i]->bond_no == 2){
	431	vtx->neigh[i]->energy = -(vtx->neigh[i]->bond[0]->xvtx->neigh[i]->bond[1]->x + vtx->neigh[i]->bond[0]->yvtx->neigh[i]->bond[1]->y + vtx->neigh[i]->bond[0]->z*vtx->neigh[i]->bond[1]->z)/vtx->neigh[i]->bond[0]->bond_length/vtx->neigh[i]->bond[1]->bond_length;
	432	delta_energy += vtx->neigh[i]->energy - backupneigh[i].energy;
	433	}
58230a	434	}
M	435
b30f45	436	// poly->k is filament persistence length (in units l_min)
M	437	delta_energy *= poly->k;
58230a	438
M	439	if(delta_energy>=0){
	440	#ifdef TS_DOUBLE_DOUBLE
	441	if(exp(-delta_energy)< drand48() )
	442	#endif
	443	#ifdef TS_DOUBLE_FLOAT
	444	if(expf(-delta_energy)< (ts_float)drand48())
	445	#endif
	446	#ifdef TS_DOUBLE_LONGDOUBLE
	447	if(expl(-delta_energy)< (ts_ldouble)drand48())
	448	#endif
	449	{
	450	//not accepted, reverting changes
	451	vtx=memcpy((void )vtx,(void )&backupvtx,sizeof(ts_vertex));
b30f45	452	for(i=0;i<vtx->neigh_no;i++){
M	453	memcpy(vtx->neigh[i],&backupneigh[i],sizeof(ts_vertex));
	454	}
58230a	455	for(i=0;i<vtx->bond_no;i++){
b30f45	456	vtx->bond[i]=memcpy((void )vtx->bond[i],(void )&backupbond[i],sizeof(ts_bond));
58230a	457	}
M	458
	459	return TS_FAIL;
	460	}
	461	}
	462
b30f45	463
58230a	464	// oldcellidx=vertex_self_avoidance(vesicle, &backupvtx[0]);
M	465	if(vtx->cell!=vesicle->clist->cell[cellidx]){
	466	retval=cell_add_vertex(vesicle->clist->cell[cellidx],vtx);
	467	// if(retval==TS_SUCCESS) cell_remove_vertex(vesicle->clist->cell[oldcellidx],vtx);
	468	if(retval==TS_SUCCESS) cell_remove_vertex(backupvtx.cell,vtx);
	469	}
	470	// if(oldcellidx);
	471	//END MONTE CARLOOOOOOO
	472	return TS_SUCCESS;
	473	}