- Posts: 4
Bug report for elastic networks in martinize.py
- buuhuu
- Topic Author
- Offline
- Fresh Boarder
Less
More
8 years 8 months ago #4852
by buuhuu
Bug report for elastic networks in martinize.py was created by buuhuu
Hi,
i noticed a small bug in the martinize.py when it comes to applying elastic networks. When assigning the the -el flag for lower cutoff length of the rubberbands, the script does not correctly assign them in the itp file.
IMHO the flawed part of code looks like:
Adding a condition for a the lower cutoff worked for my problem, but i can not assure that it will always work. Tsjerk might have a closer look at it. Here is a quick fix Peter and i worked out yesterday during the MARTINI workshop:
Let me also take this chance to thank the organizers for this beautiful workshop. It is a great experience!
Eric
i noticed a small bug in the martinize.py when it comes to applying elastic networks. When assigning the the -el flag for lower cutoff length of the rubberbands, the script does not correctly assign them in the itp file.
IMHO the flawed part of code looks like:
def rubberBands(atomList,lowerBound,upperBound,decayFactor,decayPower,forceConstant,minimumForce):
out = []
u2 = upperBound**2
while len(atomList) > 3:
bi,xi = atomList.pop(0)
for bj,xj in atomList[2:]:
# Mind the nm/A conversion -- This has to be standardized! Global use of nm?
d2 = distance2(xi,xj)/100
if d2 < u2:
dij = math.sqrt(d2)
fscl = decayFunction(dij,lowerBound,decayFactor,decayPower)
if fscl*forceConstant > minimumForce:
out.append({"atoms":(bi,bj),"parameters": (dij,"RUBBER_FC*%f"%fscl)})
return out
Adding a condition for a the lower cutoff worked for my problem, but i can not assure that it will always work. Tsjerk might have a closer look at it. Here is a quick fix Peter and i worked out yesterday during the MARTINI workshop:
def rubberBands(atomList,lowerBound,upperBound,decayFactor,decayPower,forceConstant,minimumForce):
out = []
u2 = upperBound**2
l2 = lowerBound**2
while len(atomList) > 3:
bi,xi = atomList.pop(0)
for bj,xj in atomList[2:]:
# Mind the nm/A conversion -- This has to be standardized! Global use of nm?
d2 = distance2(xi,xj)/100
if l2 < d2 < u2:
dij = math.sqrt(d2)
fscl = decayFunction(dij,lowerBound,decayFactor,decayPower)
if fscl*forceConstant > minimumForce:
out.append({"atoms":(bi,bj),"parameters": (dij,"RUBBER_FC*%f"%fscl)})
return out
Let me also take this chance to thank the organizers for this beautiful workshop. It is a great experience!
Eric
Please Log in or Create an account to join the conversation.
Time to create page: 0.086 seconds