search

NoskovLab / BROMOC

Polymer Dynamics in nanopores.
GNU General Public License v3.0
4 stars 0 forks source link

Problem running bromoc #1

Open krlitros87 opened 7 years ago

krlitros87 commented 7 years ago

Dear all, I'm a beginner in Brownian dynamics (have some experiences in molecular dynamic simulations though), and I had some problems just creating the system/files in order to use bromoc. Based on the previous I found really helpful the charmm-gui module to prepare the BD systems-> http://www.charmm-gui.org/?doc=input/gcmcbd (which of course I used to create the input files for the simulation). Surprisingly, pb-pnp worked well and I was able to create all 3 (phiv, phix and rfpar) .bin files that will be use in the gcmc/bd simulations. PHIV

  • GENERATED BY CHARMM-GUI (http://www.charmm-gui.org) v1.8 on Apr, 24. 2017. JOBID=1493041364
  • Calculate ion accessible region: phiv

OPEN read unit 10 name step3_stern_radii.crd OPEN read unit 11 name step3_pb_charge.crd PROTEIN_STRUCTURE charmm runit 10 cunit 11

ITERATION maxpnp 100000 maxphi 100000 maxcion 100000 tolphi 2.0D-6 - tolcion 5.0D-11 lambda 0.3

PHYSICAL_PARAMETERS epsp 2 epsw 80.0 temp 300.0 conc 0.15 watr 1.4 ionr 0.0

! user input MEMBRANE tmemb 30 epsm 2 htmemb 0.0 epsh 2.0 zmemb 0.0 - vmemb 0

! user input OBJECTS rsphe 0.0 xsphe 0.0 ysphe 0.0 zsphe 0.0 epss 0.0 - ax 29.016178 by 29.016178 hecyln 30 xecyln 0.0 yecyln 0.0 zecyln 0.0 epsec 80.0 eckap - bxmax 0.0 bymax 0.0 bzmax 0.0 - bxmin 0.0 bymin 0.0 bzmin 0.0 epsb 0.0 bkappa

GRID nclx 271 ncly 239 nclz 281 dcel 0.5 - xbcen 0.0 ybcen 0.0 zbcen 6.52

PBEQ zerobp reentrance COUNTERION

OPEN write file unit 12 name step4_phiv.bin WRITE kappa2 unit 12

STOP

PHIX

GENERATED BY CHARMM-GUI (http://www.charmm-gui.org) v1.8 on Apr, 24. 2017. JOBID=1493041364

  • Calculate electrostatic potential map

OPEN read unit 10 name step3_pb_radii.crd OPEN read unit 11 name step3_pb_charge.crd PROTEIN_STRUCTURE charmm runit 10 cunit 11

ITERATION maxpnp 100000 maxphi 100000 maxcion 100000 tolphi 2.0D-6 - tolcion 5.0D-11 lambda 0.2

PHYSICAL_PARAMETERS epsp 2 epsw 80.0 temp 300.0 conc 0.15 watr 0.0 ionr 0.0

MEMBRANE tmemb 30 epsm 2 htmemb 0.0 epsh 2.0 zmemb 0.0 - vmemb 0

OBJECTS rsphe 0.0 xsphe 0.0 ysphe 0.0 zsphe 0.0 epss 0.0 - ax 29.016178 by 29.016178 hecyln 30 xecyln 0.0 yecyln 0.0 zecyln 0.0 epsec 80.0 - bxmax 67.5 bymax 59.5 bzmax 76.52 - bxmin -67.5 bymin -59.5 bzmin -63.48 epsb 0.0

GRID nclx 271 ncly 239 nclz 281 dcel 1.5 - xbcen 0.0 ybcen 0.0 zbcen 6.52

PBEQ PBEQ phi nonlinear underrelax

GRID nclx 271 ncly 239 nclz 281 dcel 0.5 - xbcen 0.0 ybcen 0.0 zbcen 6.52

PBEQ phifocus pnclx 271 pncly 239 pnclz 281 pdcel 1.5 - pxbcen 0.0 pybcen 0.0 pzbcen 6.52

PBEQ phi nonlinear underrelax WRITE phi card unit 6 zfirst -200.0 zlast 200.0 COUNTERION nonlinear

OPEN write file unit 12 name step4_phix.bin WRITE phi file unit 12

STOP

RFPAR

  • GENERATED BY CHARMM-GUI (http://www.charmm-gui.org) v1.8 on Apr, 24. 2017. JOBID=1493041364
  • Calculate reaction field parameter

IONTYPE POT charge 1.0 radius 0.5 END

OPEN read unit 10 name step3_pb_radii.crd OPEN read unit 11 name step3_pb_charge.crd PROTEIN_STRUCTURE CHARMM runit 10 cunit 11

PHYSICAL_PARAMETERS epsp 2 epsw 80.0 temp 300.0 conc 0.15 watr 0.0 ionr 0.0

MEMBRANE tmemb 30 epsm 2 htmemb 0.0 epsh 2.0 zmemb 0.0 - vmemb 0

OBJECTS rsphe 0.0 xsphe 0.0 ysphe 0.0 zsphe 0.0 epss 0.0 - ax 29.016178 by 29.016178 hecyln 30 xecyln 0.0 yecyln 0.0 zecyln 0.0 epsec 80.0 - bxmax 67.5 bymax 59.5 bzmax 76.52 - bxmin -67.5 bymin -59.5 bzmin -63.48 epsb 0.0

GRID nclx 271 ncly 239 nclz 281 dcel 0.5 - xbcen 0.0 ybcen 0.0 zbcen 6.52

RFPAR box

OPEN write unit 15 file name step4_rfpar.bin WRITE rfpar POT unit 15

STOP

Unfortunately, when I used bromoc I found some issues: GCMC/BD

  • GENERATED BY CHARMM-GUI (http://www.charmm-gui.org) v1.8 on Apr, 24. 2017. JOBID=1493041364
  • GCMC/BD input file without reaction field

PARTICLES POT eps 0.087 sigma 3.142645 charge 1.0 diffusion 0.196 CLA eps 0.150 sigma 4.044680 charge -1.0 diffusion 0.203 END

!! Excess chemical potential mu should be determined BUFFERS LX 135 LY 119 LZ 140 Zcenter 6.52 POT conc 0.15 mu -0.21665465376954 voltage 0.000 LZmin -63.48 LZmax -58.48 bufferbias kb 1.0 CLA conc 0.15 mu -0.21956212105006 voltage 0.000 LZmin -63.48 LZmax -58.48 bufferbias kb 1.0 POT conc 0.15 mu -0.21665465376954 voltage 0.000 LZmin 71.52 LZmax 76.52 bufferbias kb 1.0 CLA conc 0.15 mu -0.21956212105006 voltage 0.000 LZmin 71.52 LZmax 76.52 bufferbias kb 1.0 END

!! SRPMF: Ion type dependent SRPMF POT POT c0 -0.60 c1 4.40 c2 0.90 c3 0.80 c4 0.25 CLA CLA c0 -0.50 c1 4.90 c2 0.90 c3 0.80 c4 0.25 POT CLA c0 -3.70 c1 2.90 c2 0.90 c3 0.75 c4 0.00 END

!! If model not defined, use bulk value

! SWTICHING region DIFFUSE porelength 30 pcenter 0 POT lowfraction 0.5 switchlength 10 unit 13 CLA lowfraction 0.5 switchlength 10 unit 14 END

OPEN read unit 11 file name step4_phiv.bin OPEN read unit 31 file name step4_phix.bin OPEN read unit 32 file name step4_rfpar.bin ! unit for rfpar should increase by one for the next ion species ! use same rfpar for POT and CLA OPEN read unit 33 file name step4_rfpar.bin

GSBP svdw 50 accessunit 11 phixunit 31 rfparunit 32 vzmin -63.48 vzmax 76.52

SIMULATION ncycle 10000000 ngrand 1 iseed 1493041364 nsteps 1 dt 0.01 zcont 0.0 - temp 300 nprint 100000 rho

STOP

First of, I noticed that in bromoc one should add a SYSTEM flag. Since charmm-gui didn't add it, I added only basic information of the system (SYSTEM LX 135 LY 119 LZ 126.96). After adding this line I got the following error message:

OPEN read unit 31 file name step4_phix.bin

fatal error at shell_simul : unit is greater than maxopen

STOP 1 If I'm not wrong, the unit flag only determines the ''order'' to write files, right? (please let me know if I'm wrong). Anyway, based on the previously described, I change the unit values of the files:

OPEN read unit 11 file name step4_phiv.bin OPEN read unit 12 file name step4_phix.bin OPEN read unit 13 file name step4_rfpar.bin ! unit for rfpar should increase by one for the next ion species ! use same rfpar for POT and CLA OPEN read unit 14 file name step4_rfpar.bin

But now I found the following error:

OPEN read unit 11 file name step4_phiv.bin

  open file step4_phiv.bin as unit   2

OPEN read unit 12 file name step4_phix.bin

  open file step4_phix.bin as unit   3

OPEN read unit 13 file name step4_rfpar.bin

  open file step4_rfpar.bin as unit   4

OPEN read unit 14 file name step4_rfpar.bin

  open file step4_rfpar.bin as unit   7

At line 234 of file shell_simul.f90 (unit = 7) Fortran runtime error: File already opened in another unit I know what the error means, but in the input generated by charmm-gui is stipulate that I should the same rfpar for both ions (POT and CLA), so I really don't know what to do. Probably my problem is really silly, but this is my first time working with any related to BD simulations and even with fortran, so any advices are more than welcome. Best regards, Carlos

krlitros87 commented 7 years ago

???

pablodebiase commented 7 years ago

BROMOC is a fork of GCMC/BD but it is not the same program. BROMOC has its own manual it would be worth reading it. There is a tutorial also that will help you generate new working input files. The last version of bromoc is not updated in this page. You should go to my github site. https://github.com/pablodebiase/bromoc Here is the tutorial I was talking about: https://github.com/pablodebiase/bromoc-tutorial Here is the manual: https://github.com/pablodebiase/bromoc/blob/master/cg-gcmc-bd/doc/bromoc.doc

If you want to simulate other molecules (coarse grain or all atoms) instead of just a DNA you can use the new fork called BROMOC-E. You can model any molecule as a solute (BD) or solvent (GCMC/BD) using a CHARMM force field file format. Here is the link to the code: https://github.com/pablodebiase/bromoc-e_suite

pablodebiase commented 7 years ago

Answering your previous question the unit flag is just the number assigned to the opened file. It has nothing to do with the order. The reason why it wasn't working for you was because your syntax for loading rfpar maps is incorrect. If you go to the manual:

[GSBP] {adjust} {phiv} svdw (real) vzmin (real) vzmax (real) {bspline} {nmcden} phivunit (integer) [thold27 (real)] [thold8 (real)] [repwalls (integer)] {phix} phixunit (integer) {rfpar} rfparunit (integer,integer,...) [rfparfac (real)] {rfpsingle}

.....

Here is an example of how you should load your reaction field parameters:

example 2: same repulsive potentials [phiv] for the ion types, static field [phix], and different reaction field parameters for the ion types

open read unit 11 file name phiv.bin ! particle open read unit 13 file name phix.bin open read unit 14 file name rfpar_pot.bin ! particle POT should be defined first open read unit 15 file name rfpar_cla.bin ! particle CLA should be defined second

gsbp svdw 50.0 phiv phivunit 11 phix phixunit 13 rfpar rfparunit 14,15

close unit 11 close unit 13 close unit 14 close unit 15