from cpolymer.constrain import Box,Point
from cpolymer.polymer import Polymer
from cpolymer.lsimu import LSimu
Sim = LSimu()
box = Box([-50,-50,-100],[50,50,100])
Sim.add_box(box)
N=35
step=5
start = Point(index=0,position=[-step,-step,-5])
end = Point(index=N-1,position=[0,0,12])
angle_bond=True
virtual_lp=4
P1 = Polymer(N=N,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True,virtual_lp=virtual_lp)
start.position=[-step,step,-5]
P2 = Polymer(N=N,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True,virtual_lp=virtual_lp)
start.position=[step,step,-5]
P3 = Polymer(N=N,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True,virtual_lp=virtual_lp)
start.position=[step,-step,-5]
P4 = Polymer(N=N,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True,virtual_lp=virtual_lp)
Sim.add(P1)
Sim.add(P2)
Sim.add(P3)
Sim.add(P4)
Na=15
for i in range(Na):
Sim.add_extra_bond(mol1=[1,N-i],mol2=[4,N-i],typeb=1)
for mol in range(1,4):
Sim.add_extra_bond(mol1=[mol,N-i],mol2=[mol+1,N-i],typeb=1)#molecules[-1].extrabond.append([i,1,N-i,4*N-i])
Nld=9
start = Point(index=0,position=[-0.5*step,-0.5*step,-3])
end = Point(index=Nld-1,position=[-0.5*step,0.5*step,-3])
L = Polymer(N=Nld,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True)
Sim.add(L)
start = Point(index=0,position=[-0.5*step,0.5*step,-3])
end = Point(index=Nld-1,position=[0.5*step,0.5*step,-3])
L = Polymer(N=Nld,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True)
Sim.add(L)
start = Point(index=0,position=[0.5*step,0.5*step,-3])
end = Point(index=Nld-1,position=[0.5*step,-0.5*step,-3])
L = Polymer(N=Nld,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True)
Sim.add(L)
start = Point(index=0,position=[0.5*step,-0.5*step,-3])
end = Point(index=Nld-1,position=[-0.5*step,-0.5*step,-3])
L = Polymer(N=Nld,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True)
Sim.add(L)
Sim.add_extra_bond(mol1=[5,1],mol2=[1,10],typeb=1)
Sim.add_extra_bond(mol1=[5,Nld],mol2=[2,10],typeb=1)
Sim.add_extra_bond(mol1=[6,1],mol2=[2,10],typeb=1)
Sim.add_extra_bond(mol1=[6,Nld],mol2=[3,10],typeb=1)
Sim.add_extra_bond(mol1=[7,1],mol2=[3,10],typeb=1)
Sim.add_extra_bond(mol1=[7,Nld],mol2=[4,10],typeb=1)
Sim.add_extra_bond(mol1=[8,1],mol2=[4,10],typeb=1)
Sim.add_extra_bond(mol1=[8,Nld],mol2=[1,10],typeb=1)
Nld=5
start = Point(index=0,position=[-0.2*step,-0.2*step,-1])
end = Point(index=Nld-1,position=[-0.2*step,0.2*step,-1])
L = Polymer(N=Nld,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True)
Sim.add(L)
start = Point(index=0,position=[-0.2*step,0.2*step,-1])
end = Point(index=Nld-1,position=[0.2*step,0.2*step,-1])
L = Polymer(N=Nld,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True)
Sim.add(L)
start = Point(index=0,position=[0.2*step,0.2*step,-1])
end = Point(index=Nld-1,position=[0.2*step,-0.2*step,-1])
L = Polymer(N=Nld,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True)
Sim.add(L)
start = Point(index=0,position=[0.2*step,-0.2*step,-1])
end = Point(index=Nld-1,position=[-0.2*step,-0.2*step,-1])
L = Polymer(N=Nld,lconstrain=[start,end],gconstrain=[box],angle_def={"1-1-1":[10,1]},angle_bond=True)
Sim.add(L)
Sim.add_extra_bond(mol1=[9,1],mol2=[1,15],typeb=1)
Sim.add_extra_bond(mol1=[9,Nld],mol2=[2,15],typeb=1)
Sim.add_extra_bond(mol1=[10,1],mol2=[2,15],typeb=1)
Sim.add_extra_bond(mol1=[10,Nld],mol2=[3,15],typeb=1)
Sim.add_extra_bond(mol1=[11,1],mol2=[3,15],typeb=1)
Sim.add_extra_bond(mol1=[11,Nld],mol2=[4,15],typeb=1)
Sim.add_extra_bond(mol1=[12,1],mol2=[4,15],typeb=1)
Sim.add_extra_bond(mol1=[12,Nld],mol2=[1,15],typeb=1)
rep = "./"
Sim.add_bond(typeb="harmonic",idbond=1,K=60,R0=1.)
Sim.add_pair(typep="lj/cut",idpair1=1,idpair2=1,epsilon=1,sigma=10,cutoff1=1.12)
Sim.add_angle(typea="harmonic",idangle=1,K=10,theta=180)
Sim.generate_xyz(rep + "mix.xyz",Mass="one")
Sim.generate_interactions(rep + "interactions")
Sim.generate_pdb(rep + "mix.pdb")
Sim.generate_script(rep+"basic.txt",run_length=5000,samplingrate=10,initconf="mix.xyz",
outtraj="out.dcd",outfile="out.xyz",interactions="interactions",particle="1 2"
)
#Sim.run(script="basic.txt")