Open kurokawaikki opened 2 years ago
CharlyEmpereurmot
commented
2 years ago Hi,
This is due to an update in MDAnalysis, the newer versions do NOT have this pbc argument anymore for the function center(), that was always set to None. So you can just remove this argument in the code and it should work fine. By this I mean, remove all the occurences of pbc=None. I will push an updated version of the package soon.
kurokawaikki
commented
2 years ago Thank you very much for your prompt reply. I remove all the pbc=None in the code and it worked. But I am wonder how to make a proper initial itp file. I encounter this error afterwards... ValueError: AtomGroup was not contiguous from bonds, process failed
The itp file that I used is shown in the following. Since my polymer is not symmetric, maybe I don't need the group. And my polymer is PEG-Leu. The initial values were derived with bartender. Could you give me some advice on this? Thank you very much.
[ moleculetype ]
; Name nrexcl
PGLE 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 SN3r 0 PGLE C1 1 0
2 SN3r 0 PGLE C2 2 0
3 SN3r 0 PGLE C3 3 0
4 SN3r 0 PGLE C4 4 0
5 SN3r 0 PGLE C5 5 0
6 SN2r 0 PGLE C6 6 0
7 P2 0 PGLE L61 7 0
8 TN6d 0 PGLE L62 8 1
9 C1 0 PGLE L63 9 0
10 SP2 0 PGLE C7 10 0
11 P2 0 PGLE L71 11 0
12 TN6d 0 PGLE L72 12 1
13 C1 0 PGLE L73 13 0
[bonds] ; i j funct length force.c. 1 2 1 0.35612 2072.45
2 3 1 0.297 1965.26
3 4 1 0.304 2010.50
4 5 1 0.310 1672.18
5 6 1 0.336 1370.19
6 7 1 0.309 2139.18
7 8 1 0.346 15362.05
8 9 1 0.350 6271.65
6 10 1 0.263 3604.95
10 11 1 0.271 4276.60
11 12 1 0.344 14041.80
[angles] ; i j k funct angle force_constant 1 2 3 1 130.30 16.94
2 3 4 1 128.54 16.93
3 4 5 1 127.16 19.55
4 5 6 1 123.17 16.89
5 6 7 1 70.6237 85.25
6 7 8 1 127.386 21.75
7 8 9 1 88.1504 60.27
5 6 10 1 126.06 16.37
7 6 10 1 122.15 13.52
6 10 11 1 85.1168 80.93
10 11 12 1 126.83 35.84
11 12 13 1 80.9223 76.02
CharlyEmpereurmot
commented
2 years ago The groups can be used for reducing the number of parameters to be optimized, while at the same time enhancing the consistency and physical relevance of the model. For example, in between the beads 1-6 all the bonds should in principle use the same equilibrium value and force constant -- because they represent the same chemical entities + the molecular context is quite similar -- so you can group them and the software will understand these bonded elements are related / identical parameters. You can see already that the parameters of your current model are pretty close for the bonded parameters of those bonds, which confirms this idea. Same principle applies to angles in between beads 1-6.
SwarmCG, however, is not designed for helping us to decide about how bonded elements (bonds/angles/...) should be grouped or not. How the groups are formed in your ITP is entirely your choice.
Then, about the error you encounter, this seems due to your ITP definition. I believe the bead number 13 is not bonded to the rest of the molecule (beads 1-12). The error message is not super clear about that, I will correct this!
kurokawaikki
commented
2 years ago Thank you very much!!
I solved the problem and I got the optimized parameter within 4 hours.
This is the final fitted figure:
I am wondering is this fitting good enough in your opinion. And, usually, what will be the approximate total mismatch score for a good fitted model? And will the handling of virtual site be added in the near future? Thank you very much for your great help and works.
kurokawaikki
commented
2 years ago Dear Dr. @CharlyEmpereurmot ,
I have encountered another error during the SwarmCG run.
I changed the mapping of CG for my polymer. And the optimization process keeps failing. There are two kinds of error messages during the MD and equiliberation run:
Program: gmx mdrun, version 2022 Source file: src/gromacs/gpu_utils/cudautils.cuh (line 190)
Fatal error: Unexpected cudaStreamQuery failure. CUDA error #700 (cudaErrorIllegalAddress): an illegal memory access was encountered.
2. Step 186260 Warning: pressure scaling more than 1%, mu: inf inf inf
Program: gmx mdrun, version 2022 Source file: src/gromacs/mdlib/sim_util.cpp (line 554)
Fatal error: Step 186648: The total potential energy is nan, which is not finite. The LJ and electrostatic contributions to the energy are 0 and 0, respectively. A non-finite potential energy can be caused by overlapping interactions in bonded interactions or very large or Nan coordinate values. Usually this is caused by a badly- or non-equilibrated initial configuration, incorrect interactions or parameters in the topology.
How do I solve the problem? Could you give me some hints? Thank you very much.
The new CG itp file is:
[ moleculetype ]
; Name nrexcl
PGLE 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 SN3r 0 PGLE C1 1 0
2 SN3r 0 PGLE C2 2 0
3 SN3r 0 PGLE C3 3 0
4 SN3r 0 PGLE C4 4 0
5 SN3r 0 PGLE C5 5 0
6 SN2r 0 PGLE P1 6 0
7 P2 0 PGLE L11 7 0
8 Q5 0 PGLE L12 8 1 ;should change
9 SN2r 0 PGLE P2 9 0
10 P2 0 PGLE L21 10 0
11 Q5 0 PGLE L22 11 1 ;should change
12 SN2r 0 PGLE P3 12 0
13 P2 0 PGLE L31 13 0
14 Q5 0 PGLE L32 14 1 ;should change
15 SP2 0 PGLE P4 15 0
16 P2 0 PGLE L41 16 0
17 Q5 0 PGLE L42 17 1 ;should change
[bonds]
; i j funct length force.c.
1 2 1 0.310 2056.94 ; rmsd: 1.64
2 3 1 0.305 1885.29 ; rmsd: 1.24
3 4 1 0.304 1915.49 ; rmsd: 1.48
4 5 1 0.303 1676.81 ; rmsd: 1.26
5 6 1 0.291 2710.68 ; rmsd: 2.76
6 7 1 0.317 3880.87 ; rmsd: 2.21
7 8 1 0.403 2612.33 ; rmsd: 2.06
6 9 1 0.299 3274.62 ; rmsd: 3.78
9 10 1 0.314 3347.64 ; rmsd: 2.06
10 11 1 0.403 2578.16 ; rmsd: 1.94
9 12 1 0.284 3210.37 ; rmsd: 2.71
12 13 1 0.301 3793.09 ; rmsd: 1.57
13 14 1 0.403 2612.92 ; rmsd: 2.02
12 15 1 0.262 5647.93 ; rmsd: 3.13
15 16 1 0.279 4414.50 ; rmsd: 2.26
16 17 1 0.404 2631.65 ; rmsd: 1.96
[angles] ; i j k funct angle force_constant 1 2 3 1 127.67 16.85 ; rmsd: 0.90
2 3 4 1 128.39 19.36 ; rmsd: 0.85
3 4 5 1 130.41 19.18 ; rmsd: 0.89
4 5 6 1 130.46 16.81 ; rmsd: 0.97
5 6 7 1 85.84 97.07 ; rmsd: 0.61
6 7 8 1 140.50 38.37 ; rmsd: 1.33
7 6 9 1 133.17 28.57 ; rmsd: 0.95
5 6 9 1 128.43 25.26 ; rmsd: 1.21
6 9 10 1 79.44 134.86 ; rmsd: 0.97
9 10 11 1 142.20 34.07 ; rmsd: 1.64
6 9 12 1 130.56 28.50 ; rmsd: 1.62
10 9 12 1 129.16 28.45 ; rmsd: 1.29
9 12 15 1 129.13 27.95 ; rmsd: 2.30
9 12 13 1 85.69 108.26 ; rmsd: 1.20
12 13 14 1 130.52 22.99 ; rmsd: 1.37
13 12 15 1 118.98 20.22 ; rmsd: 1.29
12 15 16 1 94.71 83.58 ; rmsd: 1.40
15 16 17 1 127.93 37.08 ; rmsd: 1.05
webridging
commented
2 months ago Hi,
When I tried to run the program with my polymers, the program gives out the error "TypeError: center() got an unexpected keyword argument 'pbc'". How can I solve the problem? Thank you.
### 做如下修改好像也行吧，不知道对不对：把pbc=None改成wrap=None。下面是MDanalysis给的注释： def _pbc_to_wrap(function): """Raises deprecation warning if 'pbc' is set and assigns value to 'wrap'""" @functools.wraps(function) def wrapped(group, *args, **kwargs): if kwargs.get('pbc', None) is not None: warnings.warn("The 'pbc' kwarg has been deprecated and will be " "removed in version 3.0., " "please use 'wrap' instead", DeprecationWarning) kwargs['wrap'] = kwargs.pop('pbc')
return function(group, *args, **kwargs)
return wrapped分子轨迹处理用MDanalysis的话，好像周期性处理有问题，统计的键长有个别和盒子边长一样，比如15nm。用mdtraj处理分子轨迹更好吧？我用过mdtraj来计算键能，他的键长统计是正确的。另外pycgtool这个工具拟合平衡键长用的也是mdtraj工具处理分子轨迹，我用pycgtool来拟合MARTINI的键长也未见有误。所以用mdtraj提供的函数处理分子轨迹会更佳吧。
Hi,
When I tried to run the program with my polymers, the program gives out the error "TypeError: center() got an unexpected keyword argument 'pbc'". How can I solve the problem? Thank you.