In the first cell we are going to create two polymers with different rigidity This polymer have some constrain to satisfy:
import random
random.seed(0)
from cpolymer.create import one_polymer
from cpolymer.constrain import Box,Point
#We set three constraint:
#They are point whoose polymer must goes by:
PS = [0.2,0.2,0.2]
PM = [10,10.2,0.2]
PE = [0.2,199,19]
start = Point(index=0,position=PS)
middle = Point(index=39,position=PM)
end = Point(index=499,position=PE)
#The box is defined here:
#Be careful: if the polymer is too rigid and the constrain points
#Too close from the box walls the polymer can be hard to generate
box = Box([-200,-200,-200],[400,400,400])
coords_wlp,bonds,type_beads,ids = one_polymer(N=400,type_bead=1,liaison={"1-1":[1.0,1]},ptolerance=0,
type_polymer="linear",start_id=0,lconstrain=[start,middle,end],
gconstrain=[box],max_trial=300000,rc=2,virtual_lp=None,rigid_constrain=False)
coords_lp,bonds,type_beads,ids = one_polymer(N=500,type_bead=1,liaison={"1-1":[1.0,1]},ptolerance=0,
type_polymer="linear",start_id=0,lconstrain=[start,middle,end],
gconstrain=[box],max_trial=500000,rc=0.1,virtual_lp=8,rigid_constrain=True)
f = figure(figsize=(10,10))
for npoly,coords in enumerate([coords_wlp,coords_lp]):
ax = f.add_subplot(2,2,1+2*npoly)
t1 = (coords[1:]-coords[:-1]) / np.sqrt(np.sum((coords[1:]-coords[:-1])**2,axis=1))[:,newaxis]
ps = np.sum(t1[1:]*t1[:-1],axis=1)
angle = np.mean(np.arccos(ps) *180/3.14)
#Plot the polymer path
plot(coords[:,0],coords[:,1],"o-",label="Mean angle %.1f"%angle)
#Plot the constrain
plot(PS[0],PS[1],"o",markersize=20)
plot(PM[0],PM[1],"o",markersize=20)
plot(PE[0],PE[1],"o",markersize=20)
legend(loc="best")
ax = f.add_subplot(2,2,2+2*npoly)
bonds = np.sqrt(np.sum((coords[1:]-coords[:-1])**2,axis=1))
plot(bonds,label="Bond size %.2f"%(bonds.mean()))
ax.set_xlabel("Bond number")
ax.set_ylabel("Bond length")
legend(loc="best")
f.tight_layout()
import random
random.seed(0)
from cpolymer.create import one_polymer
from cpolymer.constrain import Box,Point
#We set three constraint:
#They are point whoose polymer must goes by:
#The box is defined here:
#Be careful: if the polymer is too rigid and the constrain points
#Too close from the box walls the polymer can be hard to generate
box = Box([-100,-100,-100],[400,400,400])
coords_wlp,bonds,type_beads,ids = one_polymer(N=400,type_bead=1,liaison={"1-1":[2.0,1]},ptolerance=0,
type_polymer="linear",start_id=0,lconstrain=[],
gconstrain=[box],max_trial=300000,rc=0.5,virtual_lp=None,rigid_constrain=True)
coords_lp,bonds,type_beads,ids = one_polymer(N=400,type_bead=1,liaison={"1-1":[1.5,1]},ptolerance=0,
type_polymer="linear",start_id=0,lconstrain=[],
gconstrain=[box],max_trial=500000,rc=1,virtual_lp=4,rigid_constrain=False,flexible_lp=True)
f = figure(figsize=(10,10))
for npoly,coords in enumerate([coords_wlp,coords_lp]):
ax = f.add_subplot(2,2,1+2*npoly)
t1 = (coords[1:]-coords[:-1]) / np.sqrt(np.sum((coords[1:]-coords[:-1])**2,axis=1))[:,newaxis]
ps = np.sum(t1[1:]*t1[:-1],axis=1)
angle = np.mean(np.arccos(ps) *180/3.14)
#Plot the polymer path
plot(coords[:,0],coords[:,1],"o-",label="Mean angle %.1f"%angle)
#Plot the constrain
legend(loc="best")
ax = f.add_subplot(2,2,2+2*npoly)
bonds = np.sqrt(np.sum((coords[1:]-coords[:-1])**2,axis=1))
plot(bonds,label="Bond size %.2f"%(bonds.mean()))
ax.set_xlabel("Bond number")
ax.set_ylabel("Bond length")
legend(loc="best")
f.tight_layout()