2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.CHAIN'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 csa if (modecalc.eq.8) then
33 csa inquire (file="fort.40",exist=file_exist)
34 csa if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
49 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
54 c print *,"Processor",myrank," leaves READRTNS"
57 C-------------------------------------------------------------------------------
58 subroutine read_control
62 implicit real*8 (a-h,o-z)
66 logical OKRandom, prng_restart
69 include 'COMMON.IOUNITS'
70 include 'COMMON.TIME1'
71 include 'COMMON.THREAD'
72 include 'COMMON.SBRIDGE'
73 include 'COMMON.CONTROL'
76 include 'COMMON.HEADER'
77 csa include 'COMMON.CSA'
78 include 'COMMON.CHAIN'
81 include 'COMMON.FFIELD'
82 include 'COMMON.SETUP'
83 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
84 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
86 character*320 controlcard
91 read (INP,'(a)') titel
92 call card_concat(controlcard)
93 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
94 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
95 call reada(controlcard,'SEED',seed,0.0D0)
96 call random_init(seed)
97 C Set up the time limit (caution! The time must be input in minutes!)
98 read_cart=index(controlcard,'READ_CART').gt.0
99 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
100 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
101 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
102 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
103 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
104 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
105 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
106 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
107 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
108 call reada(controlcard,'DRMS',drms,0.1D0)
109 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
110 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
111 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
112 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
113 write (iout,'(a,f10.1)')'DRMS = ',drms
114 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
115 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
117 call readi(controlcard,'NZ_START',nz_start,0)
118 call readi(controlcard,'NZ_END',nz_end,0)
119 call readi(controlcard,'IZ_SC',iz_sc,0)
121 safety = 60.0d0*safety
124 call reada(controlcard,"T_BATH",t_bath,300.0d0)
125 minim=(index(controlcard,'MINIMIZE').gt.0)
126 dccart=(index(controlcard,'CART').gt.0)
127 overlapsc=(index(controlcard,'OVERLAP').gt.0)
128 overlapsc=.not.overlapsc
129 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
130 searchsc=.not.searchsc
131 sideadd=(index(controlcard,'SIDEADD').gt.0)
132 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
133 outpdb=(index(controlcard,'PDBOUT').gt.0)
134 outx=(index(controlcard,'XOUT').gt.0)
135 outmol2=(index(controlcard,'MOL2OUT').gt.0)
136 pdbref=(index(controlcard,'PDBREF').gt.0)
137 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
138 indpdb=index(controlcard,'PDBSTART')
139 extconf=(index(controlcard,'EXTCONF').gt.0)
140 call readi(controlcard,'IPRINT',iprint,0)
141 call readi(controlcard,'MAXGEN',maxgen,10000)
142 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
143 call readi(controlcard,"KDIAG",kdiag,0)
144 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
145 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
146 & write (iout,*) "RESCALE_MODE",rescale_mode
147 split_ene=index(controlcard,'SPLIT_ENE').gt.0
148 if (index(controlcard,'REGULAR').gt.0.0D0) then
149 call reada(controlcard,'WEIDIS',weidis,0.1D0)
153 call reada(controlcard,"CHECKGRAD_INC",checkgrad_inc,1.0d-4)
154 if (index(controlcard,'CHECKGRAD').gt.0) then
156 if (index(controlcard,'CART').gt.0) then
158 elseif (index(controlcard,'CARINT').gt.0) then
163 elseif (index(controlcard,'THREAD').gt.0) then
165 call readi(controlcard,'THREAD',nthread,0)
166 if (nthread.gt.0) then
167 call reada(controlcard,'WEIDIS',weidis,0.1D0)
170 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
171 stop 'Error termination in Read_Control.'
173 else if (index(controlcard,'MCMA').gt.0) then
175 else if (index(controlcard,'MCEE').gt.0) then
177 else if (index(controlcard,'MULTCONF').gt.0) then
179 else if (index(controlcard,'MAP').gt.0) then
181 call readi(controlcard,'MAP',nmap,0)
182 else if (index(controlcard,'CSA').gt.0) then
183 write(*,*) "CSA not supported in this version"
186 crc else if (index(controlcard,'ZSCORE').gt.0) then
188 crc ZSCORE is rm from UNRES, modecalc=9 is available
191 cfcm else if (index(controlcard,'MCMF').gt.0) then
193 else if (index(controlcard,'SOFTREG').gt.0) then
195 else if (index(controlcard,'CHECK_BOND').gt.0) then
197 else if (index(controlcard,'TEST').gt.0) then
199 else if (index(controlcard,'MD').gt.0) then
201 else if (index(controlcard,'RE ').gt.0) then
205 lmuca=index(controlcard,'MUCA').gt.0
206 call readi(controlcard,'MUCADYN',mucadyn,0)
207 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
208 if (lmuca .and. (me.eq.king .or. .not.out1file ))
210 write (iout,*) 'MUCADYN=',mucadyn
211 write (iout,*) 'MUCASMOOTH=',muca_smooth
214 iscode=index(controlcard,'ONE_LETTER')
215 indphi=index(controlcard,'PHI')
216 indback=index(controlcard,'BACK')
217 iranconf=index(controlcard,'RAND_CONF')
218 i2ndstr=index(controlcard,'USE_SEC_PRED')
219 gradout=index(controlcard,'GRADOUT').gt.0
220 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
222 if(me.eq.king.or..not.out1file)
223 & write (iout,'(2a)') diagmeth(kdiag),
224 & ' routine used to diagonalize matrices.'
227 c--------------------------------------------------------------------------
228 subroutine read_REMDpar
232 implicit real*8 (a-h,o-z)
234 include 'COMMON.IOUNITS'
235 include 'COMMON.TIME1'
238 include 'COMMON.LANGEVIN'
240 include 'COMMON.LANGEVIN.lang0'
242 include 'COMMON.INTERACT'
243 include 'COMMON.NAMES'
245 include 'COMMON.REMD'
246 include 'COMMON.CONTROL'
247 include 'COMMON.SETUP'
249 character*320 controlcard
250 character*3200 controlcard1
251 integer iremd_m_total
253 if(me.eq.king.or..not.out1file)
254 & write (iout,*) "REMD setup"
256 call card_concat(controlcard)
257 call readi(controlcard,"NREP",nrep,3)
258 call readi(controlcard,"NSTEX",nstex,1000)
259 call reada(controlcard,"RETMIN",retmin,10.0d0)
260 call reada(controlcard,"RETMAX",retmax,1000.0d0)
261 mremdsync=(index(controlcard,'SYNC').gt.0)
262 call readi(controlcard,"NSYN",i_sync_step,100)
263 restart1file=(index(controlcard,'REST1FILE').gt.0)
264 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
265 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
266 if(max_cache_traj_use.gt.max_cache_traj)
267 & max_cache_traj_use=max_cache_traj
268 if(me.eq.king.or..not.out1file) then
269 cd if (traj1file) then
270 crc caching is in testing - NTWX is not ignored
271 cd write (iout,*) "NTWX value is ignored"
272 cd write (iout,*) " trajectory is stored to one file by master"
273 cd write (iout,*) " before exchange at NSTEX intervals"
275 write (iout,*) "NREP= ",nrep
276 write (iout,*) "NSTEX= ",nstex
277 write (iout,*) "SYNC= ",mremdsync
278 write (iout,*) "NSYN= ",i_sync_step
279 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
282 t_exchange_only=(index(controlcard,'TONLY').gt.0)
283 call readi(controlcard,"HREMD",hremd,0)
284 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
285 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
287 if(usampl.and.hremd.gt.0) then
289 & "========== ERROR: USAMPL and HREMD cannot be used together"
291 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
298 if (index(controlcard,'TLIST').gt.0) then
300 call card_concat(controlcard1)
301 read(controlcard1,*) (remd_t(i),i=1,nrep)
302 if(me.eq.king.or..not.out1file)
303 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
306 if (index(controlcard,'MLIST').gt.0) then
308 call card_concat(controlcard1)
309 read(controlcard1,*) (remd_m(i),i=1,nrep)
310 if(me.eq.king.or..not.out1file) then
311 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
314 iremd_m_total=iremd_m_total+remd_m(i)
317 write (iout,*) 'Total number of replicas ',
318 & iremd_m_total*hremd
320 write (iout,*) 'Total number of replicas ',iremd_m_total
324 if(me.eq.king.or..not.out1file)
325 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
328 c--------------------------------------------------------------------------
329 subroutine read_MDpar
333 implicit real*8 (a-h,o-z)
335 include 'COMMON.IOUNITS'
336 include 'COMMON.TIME1'
339 include 'COMMON.LANGEVIN'
341 include 'COMMON.LANGEVIN.lang0'
343 include 'COMMON.INTERACT'
344 include 'COMMON.NAMES'
346 include 'COMMON.SETUP'
347 include 'COMMON.CONTROL'
348 include 'COMMON.SPLITELE'
350 character*320 controlcard
352 call card_concat(controlcard)
353 call readi(controlcard,"NSTEP",n_timestep,1000000)
354 call readi(controlcard,"NTWE",ntwe,100)
355 call readi(controlcard,"NTWX",ntwx,1000)
356 call reada(controlcard,"DT",d_time,1.0d-1)
357 call reada(controlcard,"DVMAX",dvmax,2.0d1)
358 call reada(controlcard,"DAMAX",damax,1.0d1)
359 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
360 call readi(controlcard,"LANG",lang,0)
361 RESPA = index(controlcard,"RESPA") .gt. 0
362 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
363 ntime_split0=ntime_split
364 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
365 ntime_split0=ntime_split
366 call reada(controlcard,"R_CUT",r_cut,2.0d0)
367 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
368 rest = index(controlcard,"REST").gt.0
369 tbf = index(controlcard,"TBF").gt.0
370 call readi(controlcard,"HMC",hmc,0)
371 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
372 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
373 tnh = index(controlcard,"NOSEHOOVER96").gt.0
374 if (RESPA.and.tnh)then
375 xiresp = index(controlcard,"XIRESP").gt.0
377 call reada(controlcard,"Q_NP",Q_np,0.1d0)
378 usampl = index(controlcard,"USAMPL").gt.0
379 mdpdb = index(controlcard,"MDPDB").gt.0
380 call reada(controlcard,"T_BATH",t_bath,300.0d0)
381 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
382 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
383 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
384 if (count_reset_moment.eq.0) count_reset_moment=1000000000
385 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
386 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
387 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
388 if (count_reset_vel.eq.0) count_reset_vel=1000000000
389 large = index(controlcard,"LARGE").gt.0
390 print_compon = index(controlcard,"PRINT_COMPON").gt.0
391 rattle = index(controlcard,"RATTLE").gt.0
392 preminim = index(controlcard,"PREMINIM").gt.0
394 dccart=(index(controlcard,'CART').gt.0)
397 c if performing umbrella sampling, fragments constrained are read from the fragment file
403 if(me.eq.king.or..not.out1file) then
405 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
407 write (iout,'(a)') "The units are:"
408 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
409 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
410 & " acceleration: angstrom/(48.9 fs)**2"
411 write (iout,'(a)') "energy: kcal/mol, temperature: K"
413 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
414 write (iout,'(a60,f10.5,a)')
415 & "Initial time step of numerical integration:",d_time,
417 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
419 write (iout,'(2a,i4,a)')
420 & "A-MTS algorithm used; initial time step for fast-varying",
421 & " short-range forces split into",ntime_split," steps."
422 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
423 & r_cut," lambda",rlamb
425 write (iout,'(2a,f10.5)')
426 & "Maximum acceleration threshold to reduce the time step",
427 & "/increase split number:",damax
428 write (iout,'(2a,f10.5)')
429 & "Maximum predicted energy drift to reduce the timestep",
430 & "/increase split number:",edriftmax
431 write (iout,'(a60,f10.5)')
432 & "Maximum velocity threshold to reduce velocities:",dvmax
433 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
434 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
435 if (rattle) write (iout,'(a60)')
436 & "Rattle algorithm used to constrain the virtual bonds"
437 if (preminim .or. iranconf.gt.0) then
439 & "Initial structure will be energy-minimized"
444 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
445 call reada(controlcard,"RWAT",rwat,1.4d0)
446 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
447 surfarea=index(controlcard,"SURFAREA").gt.0
448 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
449 if(me.eq.king.or..not.out1file)then
450 write (iout,'(/a,$)') "Langevin dynamics calculation"
453 & " with direct integration of Langevin equations"
454 else if (lang.eq.2) then
455 write (iout,'(a/)') " with TINKER stochasic MD integrator"
456 else if (lang.eq.3) then
457 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
458 else if (lang.eq.4) then
459 write (iout,'(a/)') " in overdamped mode"
461 write (iout,'(//a,i5)')
462 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
465 write (iout,'(a60,f10.5)') "Temperature:",t_bath
466 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
467 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
468 write (iout,'(a60,f10.5)')
469 & "Scaling factor of the friction forces:",scal_fric
470 if (surfarea) write (iout,'(2a,i10,a)')
471 & "Friction coefficients will be scaled by solvent-accessible",
472 & " surface area every",reset_fricmat," steps."
474 c Calculate friction coefficients and bounds of stochastic forces
475 eta=6*pi*cPoise*etawat
476 if(me.eq.king.or..not.out1file)
477 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
479 gamp=scal_fric*(pstok+rwat)*eta
480 stdfp=dsqrt(2*Rb*t_bath/d_time)
482 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
483 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
485 if(me.eq.king.or..not.out1file)then
486 write (iout,'(/2a/)')
487 & "Radii of site types and friction coefficients and std's of",
488 & " stochastic forces of fully exposed sites"
489 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
491 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
492 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
496 if(me.eq.king.or..not.out1file)then
497 write (iout,'(a)') "Berendsen bath calculation"
498 write (iout,'(a60,f10.5)') "Temperature:",t_bath
499 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
501 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
502 & count_reset_moment," steps"
504 & write (iout,'(a,i10,a)')
505 & "Velocities will be reset at random every",count_reset_vel,
508 else if (tnp .or. tnp1 .or. tnh) then
509 if (tnp .or. tnp1) then
510 write (iout,'(a)') "Nose-Poincare bath calculation"
511 if (tnp) write (iout,'(a)')
512 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
513 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
515 write (iout,'(a)') "Nose-Hoover bath calculation"
516 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
526 WDTI(i) = 1.0*d_time/nresn
533 write (iout,'(a)') "NVT-XI-RESPA algorithm"
535 write (iout,'(a)') "NVT-XO-RESPA algorithm"
538 WDTIi(i) = 1.0*d_time/nresn/ntime_split
546 write (iout,'(a60,f10.5)') "Temperature:",t_bath
547 write (iout,'(a60,f10.5)') "Q =",Q_np
549 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
550 & count_reset_moment," steps"
552 & write (iout,'(a,i10,a)')
553 & "Velocities will be reset at random every",count_reset_vel,
556 else if (hmc.gt.0) then
557 write (iout,'(a)') "Hybrid Monte Carlo calculation"
558 write (iout,'(a60,f10.5)') "Temperature:",t_bath
559 write (iout,'(a60,i10)')
560 & "Number of MD steps between Metropolis tests:",hmc
563 if(me.eq.king.or..not.out1file)
564 & write (iout,'(a31)') "Microcanonical mode calculation"
566 if(me.eq.king.or..not.out1file)then
567 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
569 write(iout,*) "MD running with constraints."
570 write(iout,*) "Equilibration time ", eq_time, " mtus."
571 write(iout,*) "Constraining ", nfrag," fragments."
572 write(iout,*) "Length of each fragment, weight and q0:"
574 write (iout,*) "Set of restraints #",iset
576 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
577 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
579 write(iout,*) "constraints between ", npair, "fragments."
580 write(iout,*) "constraint pairs, weights and q0:"
582 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
583 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
585 write(iout,*) "angle constraints within ", nfrag_back,
586 & "backbone fragments."
587 write(iout,*) "fragment, weights:"
589 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
590 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
591 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
594 iset=mod(kolor,nset)+1
597 if(me.eq.king.or..not.out1file)
598 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
601 c------------------------------------------------------------------------------
604 C Read molecular data.
606 implicit real*8 (a-h,o-z)
612 include 'COMMON.IOUNITS'
615 include 'COMMON.INTERACT'
616 include 'COMMON.LOCAL'
617 include 'COMMON.NAMES'
618 include 'COMMON.CHAIN'
619 include 'COMMON.FFIELD'
620 include 'COMMON.SBRIDGE'
621 include 'COMMON.HEADER'
622 include 'COMMON.CONTROL'
623 include 'COMMON.DBASE'
624 include 'COMMON.THREAD'
625 include 'COMMON.CONTACTS'
626 include 'COMMON.TORCNSTR'
627 include 'COMMON.TIME1'
628 include 'COMMON.BOUNDS'
630 include 'COMMON.REMD'
631 include 'COMMON.SETUP'
632 character*4 sequence(maxres)
634 double precision x(maxvar)
635 character*256 pdbfile
636 character*320 weightcard
637 character*80 weightcard_t,ucase
638 dimension itype_pdb(maxres)
639 common /pizda/ itype_pdb
640 logical seq_comp,fail
641 double precision energia(0:n_ene)
648 C Read weights of the subsequent energy terms.
661 if(me.eq.king.or..not.out1file) then
662 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
663 write (iout,*) 'Current weights for processor ',
664 & me,' set ',i2set(me)
668 call card_concat(weightcard)
669 call reada(weightcard,'WLONG',wlong,1.0D0)
670 call reada(weightcard,'WSC',wsc,wlong)
671 call reada(weightcard,'WSCP',wscp,wlong)
672 call reada(weightcard,'WELEC',welec,1.0D0)
673 call reada(weightcard,'WVDWPP',wvdwpp,welec)
674 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
675 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
676 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
677 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
678 call reada(weightcard,'WTURN3',wturn3,1.0D0)
679 call reada(weightcard,'WTURN4',wturn4,1.0D0)
680 call reada(weightcard,'WTURN6',wturn6,1.0D0)
681 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
682 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
683 call reada(weightcard,'WBOND',wbond,1.0D0)
684 call reada(weightcard,'WTOR',wtor,1.0D0)
685 call reada(weightcard,'WTORD',wtor_d,1.0D0)
686 call reada(weightcard,'WANG',wang,1.0D0)
687 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
688 call reada(weightcard,'SCAL14',scal14,0.4D0)
689 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
690 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
691 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
692 call reada(weightcard,'TEMP0',temp0,300.0d0)
693 if (index(weightcard,'SOFT').gt.0) ipot=6
694 C 12/1/95 Added weight for the multi-body term WCORR
695 call reada(weightcard,'WCORRH',wcorr,1.0D0)
696 if (wcorr4.gt.0.0d0) wcorr=wcorr4
704 hweights(i,7)=wel_loc
707 hweights(i,10)=wturn6
709 hweights(i,12)=wscloc
711 hweights(i,14)=wtor_d
712 hweights(i,15)=wstrain
713 hweights(i,16)=wvdwpp
715 hweights(i,18)=scal14
716 hweights(i,21)=wsccor
721 weights(i)=hweights(i2set(me),i)
745 call card_concat(weightcard)
746 call reada(weightcard,'WLONG',wlong,1.0D0)
747 call reada(weightcard,'WSC',wsc,wlong)
748 call reada(weightcard,'WSCP',wscp,wlong)
749 call reada(weightcard,'WELEC',welec,1.0D0)
750 call reada(weightcard,'WVDWPP',wvdwpp,welec)
751 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
752 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
753 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
754 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
755 call reada(weightcard,'WTURN3',wturn3,1.0D0)
756 call reada(weightcard,'WTURN4',wturn4,1.0D0)
757 call reada(weightcard,'WTURN6',wturn6,1.0D0)
758 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
759 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
760 call reada(weightcard,'WBOND',wbond,1.0D0)
761 call reada(weightcard,'WTOR',wtor,1.0D0)
762 call reada(weightcard,'WTORD',wtor_d,1.0D0)
763 call reada(weightcard,'WANG',wang,1.0D0)
764 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
765 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
766 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
767 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
768 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
769 call reada(weightcard,'SCAL14',scal14,0.4D0)
770 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
771 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
772 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
773 call reada(weightcard,'TEMP0',temp0,300.0d0)
774 if (index(weightcard,'SOFT').gt.0) ipot=6
775 C 12/1/95 Added weight for the multi-body term WCORR
776 call reada(weightcard,'WCORRH',wcorr,1.0D0)
777 if (wcorr4.gt.0.0d0) wcorr=wcorr4
798 weights(25)=wdfa_dist
801 weights(28)=wdfa_beta
803 if(me.eq.king.or..not.out1file)
804 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
805 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
807 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
808 10 format (/'Energy-term weights (unscaled):'//
809 & 'WSCC= ',f10.6,' (SC-SC)'/
810 & 'WSCP= ',f10.6,' (SC-p)'/
811 & 'WELEC= ',f10.6,' (p-p electr)'/
812 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
813 & 'WBOND= ',f10.6,' (stretching)'/
814 & 'WANG= ',f10.6,' (bending)'/
815 & 'WSCLOC= ',f10.6,' (SC local)'/
816 & 'WTOR= ',f10.6,' (torsional)'/
817 & 'WTORD= ',f10.6,' (double torsional)'/
818 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
819 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
820 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
821 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
822 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
823 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
824 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
825 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
826 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
827 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
828 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
829 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
830 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
831 if(me.eq.king.or..not.out1file)then
832 if (wcorr4.gt.0.0d0) then
833 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
834 & 'between contact pairs of peptide groups'
835 write (iout,'(2(a,f5.3/))')
836 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
837 & 'Range of quenching the correlation terms:',2*delt_corr
838 else if (wcorr.gt.0.0d0) then
839 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
840 & 'between contact pairs of peptide groups'
842 write (iout,'(a,f8.3)')
843 & 'Scaling factor of 1,4 SC-p interactions:',scal14
844 write (iout,'(a,f8.3)')
845 & 'General scaling factor of SC-p interactions:',scalscp
847 r0_corr=cutoff_corr-delt_corr
849 aad(i,1)=scalscp*aad(i,1)
850 aad(i,2)=scalscp*aad(i,2)
851 bad(i,1)=scalscp*bad(i,1)
852 bad(i,2)=scalscp*bad(i,2)
854 call rescale_weights(t_bath)
855 if(me.eq.king.or..not.out1file)
856 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
857 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
859 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
860 22 format (/'Energy-term weights (scaled):'//
861 & 'WSCC= ',f10.6,' (SC-SC)'/
862 & 'WSCP= ',f10.6,' (SC-p)'/
863 & 'WELEC= ',f10.6,' (p-p electr)'/
864 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
865 & 'WBOND= ',f10.6,' (stretching)'/
866 & 'WANG= ',f10.6,' (bending)'/
867 & 'WSCLOC= ',f10.6,' (SC local)'/
868 & 'WTOR= ',f10.6,' (torsional)'/
869 & 'WTORD= ',f10.6,' (double torsional)'/
870 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
871 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
872 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
873 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
874 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
875 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
876 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
877 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
878 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
879 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
880 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
881 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
882 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
883 if(me.eq.king.or..not.out1file)
884 & write (iout,*) "Reference temperature for weights calculation:",
886 call reada(weightcard,"D0CM",d0cm,3.78d0)
887 call reada(weightcard,"AKCM",akcm,15.1d0)
888 call reada(weightcard,"AKTH",akth,11.0d0)
889 call reada(weightcard,"AKCT",akct,12.0d0)
890 call reada(weightcard,"V1SS",v1ss,-1.08d0)
891 call reada(weightcard,"V2SS",v2ss,7.61d0)
892 call reada(weightcard,"V3SS",v3ss,13.7d0)
893 call reada(weightcard,"EBR",ebr,-5.50D0)
894 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
896 dyn_ss_mask(i)=.false.
900 dyn_ssbond_ij(i,j)=1.0d300
903 call reada(weightcard,"HT",Ht,0.0D0)
905 ss_depth=ebr/wsc-0.25*eps(1,1)
906 Ht=Ht/wsc-0.25*eps(1,1)
907 akcm=akcm*wstrain/wsc
908 akth=akth*wstrain/wsc
909 akct=akct*wstrain/wsc
910 v1ss=v1ss*wstrain/wsc
911 v2ss=v2ss*wstrain/wsc
912 v3ss=v3ss*wstrain/wsc
914 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
917 if(me.eq.king.or..not.out1file) then
918 write (iout,*) "Parameters of the SS-bond potential:"
919 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
921 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
922 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
923 write (iout,*)" HT",Ht
924 print *,'indpdb=',indpdb,' pdbref=',pdbref
926 if (indpdb.gt.0 .or. pdbref) then
927 read(inp,'(a)') pdbfile
928 if(me.eq.king.or..not.out1file)
929 & write (iout,'(2a)') 'PDB data will be read from file ',
930 & pdbfile(:ilen(pdbfile))
931 open(ipdbin,file=pdbfile,status='old',err=33)
933 33 write (iout,'(a)') 'Error opening PDB file.'
936 c print *,'Begin reading pdb data'
945 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
946 & (crefjlee(j,i+nres),j=1,3)
949 c print *,'Finished reading pdb data'
950 if(me.eq.king.or..not.out1file)
951 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
952 & ' nstart_sup=',nstart_sup
954 itype_pdb(i)=itype(i)
958 nct=nstart_sup+nsup-1
959 call contact(.false.,ncont_ref,icont_ref,co)
962 C Following 2 lines for diagnostics; comment out if not needed
963 write (iout,*) "Before sideadd"
965 if(me.eq.king.or..not.out1file)
966 & write(iout,*)'Adding sidechains'
973 do while (fail.and.nsi.le.maxsi)
974 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
977 if(fail) write(iout,*)'Adding sidechain failed for res ',
978 & i,' after ',nsi,' trials'
981 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
984 C Following 2 lines for diagnostics; comment out if not needed
985 c write (iout,*) "After sideadd"
988 if (indpdb.eq.0) then
989 C Read sequence if not taken from the pdb file.
991 c print *,'nres=',nres
992 if (iscode.gt.0) then
993 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
995 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
997 C Convert sequence to numeric code
999 itype(i)=rescode(i,sequence(i),iscode)
1001 if (itype(2).eq.10.and.itype(1).eq.ntyp1) then
1003 & "Glycine is the first full residue, initial dummy deleted"
1009 if (itype(nres-1).eq.10.and.itype(nres).eq.ntyp1) then
1011 & "Glycine is the last full residue, terminal dummy deleted"
1015 C Assign initial virtual bond lengths
1021 vbld(i+nres)=dsc(itype(i))
1022 vbld_inv(i+nres)=dsc_inv(itype(i))
1023 c write (iout,*) "i",i," itype",itype(i),
1024 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
1028 c print '(20i4)',(itype(i),i=1,nres)
1031 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1033 if (itype(i).eq.21) then
1037 else if (itype(i+1).ne.20) then
1039 else if (itype(i).ne.20) then
1046 if(me.eq.king.or..not.out1file)then
1047 write (iout,*) "ITEL"
1049 write (iout,*) i,itype(i),itel(i)
1051 print *,'Call Read_Bridge.'
1054 C 8/13/98 Set limits to generating the dihedral angles
1059 read (inp,*) ndih_constr
1060 if (ndih_constr.gt.0) then
1062 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1063 if(me.eq.king.or..not.out1file)then
1065 & 'There are',ndih_constr,' constraints on phi angles.'
1067 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1071 phi0(i)=deg2rad*phi0(i)
1072 drange(i)=deg2rad*drange(i)
1074 if(me.eq.king.or..not.out1file)
1075 & write (iout,*) 'FTORS',ftors
1078 phibound(1,ii) = phi0(i)-drange(i)
1079 phibound(2,ii) = phi0(i)+drange(i)
1084 if (me.eq.king) then
1086 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1088 write (iout,'(a3,i5,2f10.1)')
1089 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1095 cd print *,'NNT=',NNT,' NCT=',NCT
1096 if (itype(1).eq.21) nnt=2
1097 if (itype(nres).eq.21) nct=nct-1
1099 C Bartek:READ init_vars
1100 C Initialize variables!
1101 C Juyong:READ read_info
1102 C READ fragment information!!
1103 C both routines should be in dfa.F file!!
1106 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1107 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1109 print*, 'init_dfa_vars finished!'
1111 print*, 'read_dfa_info finished!'
1116 if(me.eq.king.or..not.out1file)
1117 & write (iout,'(a,i3)') 'nsup=',nsup
1119 if (nsup.le.(nct-nnt+1)) then
1120 do i=0,nct-nnt+1-nsup
1121 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1127 & 'Error - sequences to be superposed do not match.'
1130 do i=0,nsup-(nct-nnt+1)
1131 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1133 nstart_sup=nstart_sup+i
1139 & 'Error - sequences to be superposed do not match.'
1142 if (nsup.eq.0) nsup=nct-nnt
1143 if (nstart_sup.eq.0) nstart_sup=nnt
1144 if (nstart_seq.eq.0) nstart_seq=nnt
1145 if(me.eq.king.or..not.out1file)
1146 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1147 & ' nstart_seq=',nstart_seq
1149 c--- Zscore rms -------
1150 if (nz_start.eq.0) nz_start=nnt
1151 if (nz_end.eq.0 .and. nsup.gt.0) then
1153 else if (nz_end.eq.0) then
1156 if(me.eq.king.or..not.out1file)then
1157 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1158 write (iout,*) 'IZ_SC=',iz_sc
1160 c----------------------
1163 if (.not.pdbref) then
1164 call read_angles(inp,*38)
1166 38 write (iout,'(a)') 'Error reading reference structure.'
1168 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1169 stop 'Error reading reference structure'
1173 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1182 call contact(.true.,ncont_ref,icont_ref,co)
1184 if(me.eq.king.or..not.out1file)
1185 & write (iout,*) 'Contact order:',co
1187 if(me.eq.king.or..not.out1file)
1188 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1191 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1193 if(me.eq.king.or..not.out1file)
1194 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1195 & icont_ref(1,i),' ',
1196 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1200 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1201 if (constr_dist.gt.0) then
1202 call read_dist_constr
1206 if (constr_homology.gt.0) then
1207 call read_constr_homology
1208 if (indpdb.gt.0 .or. pdbref) then
1211 c(j,i)=crefjlee(j,i)
1212 cref(j,i)=crefjlee(j,i)
1217 write (iout,*) "Array C"
1219 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1220 & (c(j,i+nres),j=1,3)
1222 write (iout,*) "Array Cref"
1224 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1225 & (cref(j,i+nres),j=1,3)
1228 call int_from_cart1(.false.)
1229 call sc_loc_geom(.false.)
1231 thetaref(i)=theta(i)
1236 dc(j,i)=c(j,i+1)-c(j,i)
1237 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1242 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1243 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1251 if (nhpb.gt.0) call hpb_partition
1252 c write (iout,*) "After read_dist_constr nhpb",nhpb
1254 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1255 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1256 & modecalc.ne.10) then
1257 C If input structure hasn't been supplied from the PDB file read or generate
1259 if (iranconf.eq.0 .and. .not. extconf) then
1260 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1261 & write (iout,'(a)') 'Initial geometry will be read in.'
1263 read (inp,*) time,potE,uconst,t_bath,
1264 & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
1265 read(inp,'(8f10.5)',end=36,err=36)
1266 & ((c(l,k),l=1,3),k=1,nres),
1267 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1268 call int_from_cart1(.false.)
1271 dc(j,i)=c(j,i+1)-c(j,i)
1272 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1276 if (itype(i).ne.10) then
1278 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1279 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1285 call read_angles(inp,*36)
1288 36 write (iout,'(a)') 'Error reading angle file.'
1290 call mpi_finalize( MPI_COMM_WORLD,IERR )
1292 stop 'Error reading angle file.'
1294 else if (extconf) then
1295 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1296 & write (iout,'(a)') 'Extended chain initial geometry.'
1298 theta(i)=90d0*deg2rad
1301 phi(i)=180d0*deg2rad
1304 alph(i)=110d0*deg2rad
1307 omeg(i)=-120d0*deg2rad
1310 if(me.eq.king.or..not.out1file)
1311 & write (iout,'(a)') 'Random-generated initial geometry.'
1315 if (me.eq.king .or. fg_rank.eq.0 .and. (
1316 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1320 call gen_rand_conf(itmp,*30)
1322 30 write (iout,*) 'Failed to generate random conformation',
1323 & ', itrial=',itrial
1324 write (*,*) 'Processor:',me,
1325 & ' Failed to generate random conformation',
1335 write (iout,'(a,i3,a)') 'Processor:',me,
1336 & ' error in generating random conformation.'
1337 write (*,'(a,i3,a)') 'Processor:',me,
1338 & ' error in generating random conformation.'
1341 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1348 elseif (modecalc.eq.4) then
1349 read (inp,'(a)') intinname
1350 open (intin,file=intinname,status='old',err=333)
1351 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1352 & write (iout,'(a)') 'intinname',intinname
1353 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1355 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1357 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1359 stop 'Error opening angle file.'
1363 C Generate distance constraints, if the PDB structure is to be regularized.
1364 if (nthread.gt.0) then
1365 call read_threadbase
1367 write (iout,*) "READRTNS: Calling setup_var"
1370 if (me.eq.king .or. .not. out1file)
1372 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1373 write (iout,'(/a,i3,a)')
1374 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1375 write (iout,'(20i4)') (iss(i),i=1,ns)
1377 write(iout,*)"Running with dynamic disulfide-bond formation"
1379 write (iout,'(/a/)') 'Pre-formed links are:'
1385 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1386 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1392 if (ns.gt.0.and.dyn_ss) then
1396 forcon(i-nss)=forcon(i)
1403 dyn_ss_mask(iss(i))=.true.
1406 if (i2ndstr.gt.0) call secstrp2dihc
1407 c call geom_to_var(nvar,x)
1408 c call etotal(energia(0))
1409 c call enerprint(energia(0))
1410 c call briefout(0,etot)
1412 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1413 cd write (iout,'(a)') 'Variable list:'
1414 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1416 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1417 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1418 & 'Processor',myrank,': end reading molecular data.'
1422 c--------------------------------------------------------------------------
1423 logical function seq_comp(itypea,itypeb,length)
1425 integer length,itypea(length),itypeb(length)
1428 if (itypea(i).ne.itypeb(i)) then
1436 c-----------------------------------------------------------------------------
1437 subroutine read_bridge
1438 C Read information about disulfide bridges.
1439 implicit real*8 (a-h,o-z)
1440 include 'DIMENSIONS'
1444 include 'COMMON.IOUNITS'
1445 include 'COMMON.GEO'
1446 include 'COMMON.VAR'
1447 include 'COMMON.INTERACT'
1448 include 'COMMON.LOCAL'
1449 include 'COMMON.NAMES'
1450 include 'COMMON.CHAIN'
1451 include 'COMMON.FFIELD'
1452 include 'COMMON.SBRIDGE'
1453 include 'COMMON.HEADER'
1454 include 'COMMON.CONTROL'
1455 include 'COMMON.DBASE'
1456 include 'COMMON.THREAD'
1457 include 'COMMON.TIME1'
1458 include 'COMMON.SETUP'
1459 C Read bridging residues.
1460 read (inp,*) ns,(iss(i),i=1,ns)
1462 if(me.eq.king.or..not.out1file)
1463 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1464 C Check whether the specified bridging residues are cystines.
1466 if (itype(iss(i)).ne.1) then
1467 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1468 & 'Do you REALLY think that the residue ',
1469 & restyp(itype(iss(i))),i,
1470 & ' can form a disulfide bridge?!!!'
1471 write (*,'(2a,i3,a)')
1472 & 'Do you REALLY think that the residue ',
1473 & restyp(itype(iss(i))),i,
1474 & ' can form a disulfide bridge?!!!'
1476 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1481 C Read preformed bridges.
1483 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1485 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1488 C Check if the residues involved in bridges are in the specified list of
1489 C bridging residues.
1492 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1493 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1494 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1495 & ' contains residues present in other pairs.'
1496 write (*,'(a,i3,a)') 'Disulfide pair',i,
1497 & ' contains residues present in other pairs.'
1499 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1505 if (ihpb(i).eq.iss(j)) goto 10
1507 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1510 if (jhpb(i).eq.iss(j)) goto 20
1512 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1518 ihpb(i)=ihpb(i)+nres
1519 jhpb(i)=jhpb(i)+nres
1525 c----------------------------------------------------------------------------
1526 subroutine read_x(kanal,*)
1527 implicit real*8 (a-h,o-z)
1528 include 'DIMENSIONS'
1529 include 'COMMON.GEO'
1530 include 'COMMON.VAR'
1531 include 'COMMON.CHAIN'
1532 include 'COMMON.IOUNITS'
1533 include 'COMMON.CONTROL'
1534 include 'COMMON.LOCAL'
1535 include 'COMMON.INTERACT'
1536 c Read coordinates from input
1538 read(kanal,'(8f10.5)',end=10,err=10)
1539 & ((c(l,k),l=1,3),k=1,nres),
1540 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1543 c(j,2*nres)=c(j,nres)
1545 call int_from_cart1(.false.)
1548 dc(j,i)=c(j,i+1)-c(j,i)
1549 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1553 if (itype(i).ne.10) then
1555 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1556 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1564 c----------------------------------------------------------------------------
1565 subroutine read_threadbase
1566 implicit real*8 (a-h,o-z)
1567 include 'DIMENSIONS'
1568 include 'COMMON.IOUNITS'
1569 include 'COMMON.GEO'
1570 include 'COMMON.VAR'
1571 include 'COMMON.INTERACT'
1572 include 'COMMON.LOCAL'
1573 include 'COMMON.NAMES'
1574 include 'COMMON.CHAIN'
1575 include 'COMMON.FFIELD'
1576 include 'COMMON.SBRIDGE'
1577 include 'COMMON.HEADER'
1578 include 'COMMON.CONTROL'
1579 include 'COMMON.DBASE'
1580 include 'COMMON.THREAD'
1581 include 'COMMON.TIME1'
1582 C Read pattern database for threading.
1583 read (icbase,*) nseq
1585 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1586 & nres_base(2,i),nres_base(3,i)
1587 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1589 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1590 c & nres_base(2,i),nres_base(3,i)
1591 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1595 if (weidis.eq.0.0D0) weidis=0.1D0
1604 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1605 write (iout,'(a,i5)') 'nexcl: ',nexcl
1606 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1609 c------------------------------------------------------------------------------
1610 subroutine setup_var
1611 implicit real*8 (a-h,o-z)
1612 include 'DIMENSIONS'
1613 include 'COMMON.IOUNITS'
1614 include 'COMMON.GEO'
1615 include 'COMMON.VAR'
1616 include 'COMMON.INTERACT'
1617 include 'COMMON.LOCAL'
1618 include 'COMMON.NAMES'
1619 include 'COMMON.CHAIN'
1620 include 'COMMON.FFIELD'
1621 include 'COMMON.SBRIDGE'
1622 include 'COMMON.HEADER'
1623 include 'COMMON.CONTROL'
1624 include 'COMMON.DBASE'
1625 include 'COMMON.THREAD'
1626 include 'COMMON.TIME1'
1627 C Set up variable list.
1633 if (itype(i).ne.10) then
1635 ialph(i,1)=nvar+nside
1639 if (indphi.gt.0) then
1641 else if (indback.gt.0) then
1646 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1649 c----------------------------------------------------------------------------
1650 subroutine gen_dist_constr
1651 C Generate CA distance constraints.
1652 implicit real*8 (a-h,o-z)
1653 include 'DIMENSIONS'
1654 include 'COMMON.IOUNITS'
1655 include 'COMMON.GEO'
1656 include 'COMMON.VAR'
1657 include 'COMMON.INTERACT'
1658 include 'COMMON.LOCAL'
1659 include 'COMMON.NAMES'
1660 include 'COMMON.CHAIN'
1661 include 'COMMON.FFIELD'
1662 include 'COMMON.SBRIDGE'
1663 include 'COMMON.HEADER'
1664 include 'COMMON.CONTROL'
1665 include 'COMMON.DBASE'
1666 include 'COMMON.THREAD'
1667 include 'COMMON.TIME1'
1668 dimension itype_pdb(maxres)
1669 common /pizda/ itype_pdb
1671 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1672 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1673 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1675 do i=nstart_sup,nstart_sup+nsup-1
1676 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1677 cd & ' seq_pdb', restyp(itype_pdb(i))
1678 do j=i+2,nstart_sup+nsup-1
1680 ihpb(nhpb)=i+nstart_seq-nstart_sup
1681 jhpb(nhpb)=j+nstart_seq-nstart_sup
1683 dhpb(nhpb)=dist(i,j)
1686 cd write (iout,'(a)') 'Distance constraints:'
1691 cd if (ii.gt.nres) then
1696 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1697 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1698 cd & dhpb(i),forcon(i)
1702 c----------------------------------------------------------------------------
1704 implicit real*8 (a-h,o-z)
1705 include 'DIMENSIONS'
1706 include 'COMMON.MAP'
1707 include 'COMMON.IOUNITS'
1708 character*3 angid(4) /'THE','PHI','ALP','OME'/
1709 character*80 mapcard,ucase
1711 read (inp,'(a)') mapcard
1712 mapcard=ucase(mapcard)
1713 if (index(mapcard,'PHI').gt.0) then
1715 else if (index(mapcard,'THE').gt.0) then
1717 else if (index(mapcard,'ALP').gt.0) then
1719 else if (index(mapcard,'OME').gt.0) then
1722 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1723 stop 'Error - illegal variable spec in MAP card.'
1725 call readi (mapcard,'RES1',res1(imap),0)
1726 call readi (mapcard,'RES2',res2(imap),0)
1727 if (res1(imap).eq.0) then
1728 res1(imap)=res2(imap)
1729 else if (res2(imap).eq.0) then
1730 res2(imap)=res1(imap)
1732 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1734 & 'Error - illegal definition of variable group in MAP.'
1735 stop 'Error - illegal definition of variable group in MAP.'
1737 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1738 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1739 call readi(mapcard,'NSTEP',nstep(imap),0)
1740 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1742 & 'Illegal boundary and/or step size specification in MAP.'
1743 stop 'Illegal boundary and/or step size specification in MAP.'
1748 c----------------------------------------------------------------------------
1749 csa subroutine csaread
1750 csa implicit real*8 (a-h,o-z)
1751 csa include 'DIMENSIONS'
1752 csa include 'COMMON.IOUNITS'
1753 csa include 'COMMON.GEO'
1754 csa include 'COMMON.CSA'
1755 csa include 'COMMON.BANK'
1756 csa include 'COMMON.CONTROL'
1757 csa character*80 ucase
1758 csa character*620 mcmcard
1759 csa call card_concat(mcmcard)
1761 csa call readi(mcmcard,'NCONF',nconf,50)
1762 csa call readi(mcmcard,'NADD',nadd,0)
1763 csa call readi(mcmcard,'JSTART',jstart,1)
1764 csa call readi(mcmcard,'JEND',jend,1)
1765 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1766 csa call readi(mcmcard,'N0',n0,1)
1767 csa call readi(mcmcard,'N1',n1,6)
1768 csa call readi(mcmcard,'N2',n2,4)
1769 csa call readi(mcmcard,'N3',n3,0)
1770 csa call readi(mcmcard,'N4',n4,0)
1771 csa call readi(mcmcard,'N5',n5,0)
1772 csa call readi(mcmcard,'N6',n6,10)
1773 csa call readi(mcmcard,'N7',n7,0)
1774 csa call readi(mcmcard,'N8',n8,0)
1775 csa call readi(mcmcard,'N9',n9,0)
1776 csa call readi(mcmcard,'N14',n14,0)
1777 csa call readi(mcmcard,'N15',n15,0)
1778 csa call readi(mcmcard,'N16',n16,0)
1779 csa call readi(mcmcard,'N17',n17,0)
1780 csa call readi(mcmcard,'N18',n18,0)
1782 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1784 csa call readi(mcmcard,'NDIFF',ndiff,2)
1785 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1786 csa call readi(mcmcard,'IS1',is1,1)
1787 csa call readi(mcmcard,'IS2',is2,8)
1788 csa call readi(mcmcard,'NRAN0',nran0,4)
1789 csa call readi(mcmcard,'NRAN1',nran1,2)
1790 csa call readi(mcmcard,'IRR',irr,1)
1791 csa call readi(mcmcard,'NSEED',nseed,20)
1792 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1793 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1794 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1795 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1796 csa call readi(mcmcard,'ICMAX',icmax,3)
1797 csa call readi(mcmcard,'IRESTART',irestart,0)
1798 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1801 csa call reada(mcmcard,'DELE',dele,20.0d0)
1802 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1803 csa call readi(mcmcard,'IREF',iref,0)
1804 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1805 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1806 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1807 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1808 csa write (iout,*) "NCONF_IN",nconf_in
1811 c----------------------------------------------------------------------------
1812 cfmc subroutine mcmfread
1813 cfmc implicit real*8 (a-h,o-z)
1814 cfmc include 'DIMENSIONS'
1815 cfmc include 'COMMON.MCMF'
1816 cfmc include 'COMMON.IOUNITS'
1817 cfmc include 'COMMON.GEO'
1818 cfmc character*80 ucase
1819 cfmc character*620 mcmcard
1820 cfmc call card_concat(mcmcard)
1822 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1823 cfmc write(iout,*)'MAXRANT=',maxrant
1824 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1825 cfmc write(iout,*)'MAXFAM=',maxfam
1826 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1827 cfmc write(iout,*)'NNET1=',nnet1
1828 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1829 cfmc write(iout,*)'NNET2=',nnet2
1830 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1831 cfmc write(iout,*)'NNET3=',nnet3
1832 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1833 cfmc write(iout,*)'ILASTT=',ilastt
1834 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1835 cfmc write(iout,*)'MAXSTR=',maxstr
1836 cfmc maxstr_f=maxstr/maxfam
1837 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1838 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1839 cfmc write(iout,*)'NMCMF=',nmcmf
1840 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1841 cfmc write(iout,*)'IFOCUS=',ifocus
1842 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1843 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1844 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1845 cfmc write(iout,*)'INTPRT=',intprt
1846 cfmc call readi(mcmcard,'IPRT',iprt,100)
1847 cfmc write(iout,*)'IPRT=',iprt
1848 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1849 cfmc write(iout,*)'IMAXTR=',imaxtr
1850 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1851 cfmc write(iout,*)'MAXEVEN=',maxeven
1852 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1853 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1854 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1855 cfmc write(iout,*)'INIMIN=',inimin
1856 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1857 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1858 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1859 cfmc write(iout,*)'NTHREAD=',nthread
1860 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1861 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1862 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1863 cfmc write(iout,*)'MAXPERT=',maxpert
1864 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1865 cfmc write(iout,*)'IRMSD=',irmsd
1866 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1867 cfmc write(iout,*)'DENEMIN=',denemin
1868 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1869 cfmc write(iout,*)'RCUT1S=',rcut1s
1870 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1871 cfmc write(iout,*)'RCUT1E=',rcut1e
1872 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1873 cfmc write(iout,*)'RCUT2S=',rcut2s
1874 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1875 cfmc write(iout,*)'RCUT2E=',rcut2e
1876 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1877 cfmc write(iout,*)'DPERT1=',d_pert1
1878 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1879 cfmc write(iout,*)'DPERT1A=',d_pert1a
1880 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1881 cfmc write(iout,*)'DPERT2=',d_pert2
1882 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1883 cfmc write(iout,*)'DPERT2A=',d_pert2a
1884 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1885 cfmc write(iout,*)'DPERT2B=',d_pert2b
1886 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1887 cfmc write(iout,*)'DPERT2C=',d_pert2c
1888 cfmc d_pert1=deg2rad*d_pert1
1889 cfmc d_pert1a=deg2rad*d_pert1a
1890 cfmc d_pert2=deg2rad*d_pert2
1891 cfmc d_pert2a=deg2rad*d_pert2a
1892 cfmc d_pert2b=deg2rad*d_pert2b
1893 cfmc d_pert2c=deg2rad*d_pert2c
1894 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1895 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1896 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1897 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1898 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1899 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1900 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1901 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1902 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1903 cfmc write(iout,*)'RCUTINI=',rcutini
1904 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1905 cfmc write(iout,*)'GRAT=',grat
1906 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1907 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1911 c----------------------------------------------------------------------------
1913 implicit real*8 (a-h,o-z)
1914 include 'DIMENSIONS'
1915 include 'COMMON.MCM'
1916 include 'COMMON.MCE'
1917 include 'COMMON.IOUNITS'
1919 character*320 mcmcard
1920 call card_concat(mcmcard)
1921 call readi(mcmcard,'MAXACC',maxacc,100)
1922 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1923 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1924 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1925 call readi(mcmcard,'MAXREPM',maxrepm,200)
1926 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1927 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1928 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1929 call reada(mcmcard,'E_UP',e_up,5.0D0)
1930 call reada(mcmcard,'DELTE',delte,0.1D0)
1931 call readi(mcmcard,'NSWEEP',nsweep,5)
1932 call readi(mcmcard,'NSTEPH',nsteph,0)
1933 call readi(mcmcard,'NSTEPC',nstepc,0)
1934 call reada(mcmcard,'TMIN',tmin,298.0D0)
1935 call reada(mcmcard,'TMAX',tmax,298.0D0)
1936 call readi(mcmcard,'NWINDOW',nwindow,0)
1937 call readi(mcmcard,'PRINT_MC',print_mc,0)
1938 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1939 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1940 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1941 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1942 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1943 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1944 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1945 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1946 if (nwindow.gt.0) then
1947 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1949 winlen(i)=winend(i)-winstart(i)+1
1952 if (tmax.lt.tmin) tmax=tmin
1953 if (tmax.eq.tmin) then
1957 if (nstepc.gt.0 .and. nsteph.gt.0) then
1958 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1959 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1961 C Probabilities of different move types
1962 sumpro_type(0)=0.0D0
1963 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1964 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1965 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1966 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1967 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1968 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1969 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1971 print *,'i',i,' sumprotype',sumpro_type(i)
1972 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1973 print *,'i',i,' sumprotype',sumpro_type(i)
1977 c----------------------------------------------------------------------------
1978 subroutine read_minim
1979 implicit real*8 (a-h,o-z)
1980 include 'DIMENSIONS'
1981 include 'COMMON.MINIM'
1982 include 'COMMON.IOUNITS'
1984 character*320 minimcard
1985 call card_concat(minimcard)
1986 call readi(minimcard,'MAXMIN',maxmin,2000)
1987 call readi(minimcard,'MAXFUN',maxfun,5000)
1988 call readi(minimcard,'MINMIN',minmin,maxmin)
1989 call readi(minimcard,'MINFUN',minfun,maxmin)
1990 call reada(minimcard,'TOLF',tolf,1.0D-2)
1991 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1992 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1993 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1994 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1995 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1996 & 'Options in energy minimization:'
1997 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1998 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1999 & 'MinMin:',MinMin,' MinFun:',MinFun,
2000 & ' TolF:',TolF,' RTolF:',RTolF
2003 c----------------------------------------------------------------------------
2004 subroutine read_angles(kanal,*)
2005 implicit real*8 (a-h,o-z)
2006 include 'DIMENSIONS'
2007 include 'COMMON.GEO'
2008 include 'COMMON.VAR'
2009 include 'COMMON.CHAIN'
2010 include 'COMMON.IOUNITS'
2011 include 'COMMON.CONTROL'
2012 c Read angles from input
2014 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
2015 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
2016 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2017 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2020 c 9/7/01 avoid 180 deg valence angle
2021 if (theta(i).gt.179.99d0) theta(i)=179.99d0
2023 theta(i)=deg2rad*theta(i)
2024 phi(i)=deg2rad*phi(i)
2025 alph(i)=deg2rad*alph(i)
2026 omeg(i)=deg2rad*omeg(i)
2031 c----------------------------------------------------------------------------
2032 subroutine reada(rekord,lancuch,wartosc,default)
2034 character*(*) rekord,lancuch
2035 double precision wartosc,default
2038 iread=index(rekord,lancuch)
2039 if (iread.eq.0) then
2043 iread=iread+ilen(lancuch)+1
2044 read (rekord(iread:),*,err=10,end=10) wartosc
2049 c----------------------------------------------------------------------------
2050 subroutine readi(rekord,lancuch,wartosc,default)
2052 character*(*) rekord,lancuch
2053 integer wartosc,default
2056 iread=index(rekord,lancuch)
2057 if (iread.eq.0) then
2061 iread=iread+ilen(lancuch)+1
2062 read (rekord(iread:),*,err=10,end=10) wartosc
2067 c----------------------------------------------------------------------------
2068 subroutine multreadi(rekord,lancuch,tablica,dim,default)
2071 integer tablica(dim),default
2072 character*(*) rekord,lancuch
2079 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2080 if (iread.eq.0) return
2081 iread=iread+ilen(lancuch)+1
2082 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2085 c----------------------------------------------------------------------------
2086 subroutine multreada(rekord,lancuch,tablica,dim,default)
2089 double precision tablica(dim),default
2090 character*(*) rekord,lancuch
2097 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2098 if (iread.eq.0) return
2099 iread=iread+ilen(lancuch)+1
2100 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2103 c----------------------------------------------------------------------------
2104 subroutine openunits
2105 implicit real*8 (a-h,o-z)
2106 include 'DIMENSIONS'
2109 character*16 form,nodename
2112 include 'COMMON.SETUP'
2113 include 'COMMON.IOUNITS'
2115 include 'COMMON.CONTROL'
2116 integer lenpre,lenpot,ilen,lentmp
2118 character*3 out1file_text,ucase
2121 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2122 call getenv_loc("PREFIX",prefix)
2124 call getenv_loc("POT",pot)
2125 call getenv_loc("DIRTMP",tmpdir)
2126 call getenv_loc("CURDIR",curdir)
2127 call getenv_loc("OUT1FILE",out1file_text)
2128 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2129 out1file_text=ucase(out1file_text)
2130 if (out1file_text(1:1).eq."Y") then
2133 out1file=fg_rank.gt.0
2138 if (lentmp.gt.0) then
2139 write (*,'(80(1h!))')
2140 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2141 write (*,'(80(1h!))')
2142 write (*,*)"All output files will be on node /tmp directory."
2144 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2145 if (me.eq.king) then
2146 write (*,*) "The master node is ",nodename
2147 else if (fg_rank.eq.0) then
2148 write (*,*) "I am the CG slave node ",nodename
2150 write (*,*) "I am the FG slave node ",nodename
2153 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2154 lenpre = lentmp+lenpre+1
2156 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2157 C Get the names and open the input files
2158 #if defined(WINIFL) || defined(WINPGI)
2159 open(1,file=pref_orig(:ilen(pref_orig))//
2160 & '.inp',status='old',readonly,shared)
2161 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2162 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2163 C Get parameter filenames and open the parameter files.
2164 call getenv_loc('BONDPAR',bondname)
2165 open (ibond,file=bondname,status='old',readonly,shared)
2166 call getenv_loc('THETPAR',thetname)
2167 open (ithep,file=thetname,status='old',readonly,shared)
2169 call getenv_loc('THETPARPDB',thetname_pdb)
2170 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2172 call getenv_loc('ROTPAR',rotname)
2173 open (irotam,file=rotname,status='old',readonly,shared)
2175 call getenv_loc('ROTPARPDB',rotname_pdb)
2176 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2178 call getenv_loc('TORPAR',torname)
2179 open (itorp,file=torname,status='old',readonly,shared)
2180 call getenv_loc('TORDPAR',tordname)
2181 open (itordp,file=tordname,status='old',readonly,shared)
2182 call getenv_loc('FOURIER',fouriername)
2183 open (ifourier,file=fouriername,status='old',readonly,shared)
2184 call getenv_loc('ELEPAR',elename)
2185 open (ielep,file=elename,status='old',readonly,shared)
2186 call getenv_loc('SIDEPAR',sidename)
2187 open (isidep,file=sidename,status='old',readonly,shared)
2188 #elif (defined CRAY) || (defined AIX)
2189 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2191 c print *,"Processor",myrank," opened file 1"
2192 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2193 c print *,"Processor",myrank," opened file 9"
2194 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2195 C Get parameter filenames and open the parameter files.
2196 call getenv_loc('BONDPAR',bondname)
2197 open (ibond,file=bondname,status='old',action='read')
2198 c print *,"Processor",myrank," opened file IBOND"
2199 call getenv_loc('THETPAR',thetname)
2200 open (ithep,file=thetname,status='old',action='read')
2201 c print *,"Processor",myrank," opened file ITHEP"
2203 call getenv_loc('THETPARPDB',thetname_pdb)
2204 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2206 call getenv_loc('ROTPAR',rotname)
2207 open (irotam,file=rotname,status='old',action='read')
2208 c print *,"Processor",myrank," opened file IROTAM"
2210 call getenv_loc('ROTPARPDB',rotname_pdb)
2211 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2213 call getenv_loc('TORPAR',torname)
2214 open (itorp,file=torname,status='old',action='read')
2215 c print *,"Processor",myrank," opened file ITORP"
2216 call getenv_loc('TORDPAR',tordname)
2217 open (itordp,file=tordname,status='old',action='read')
2218 c print *,"Processor",myrank," opened file ITORDP"
2219 call getenv_loc('SCCORPAR',sccorname)
2220 open (isccor,file=sccorname,status='old',action='read')
2221 c print *,"Processor",myrank," opened file ISCCOR"
2222 call getenv_loc('FOURIER',fouriername)
2223 open (ifourier,file=fouriername,status='old',action='read')
2224 c print *,"Processor",myrank," opened file IFOURIER"
2225 call getenv_loc('ELEPAR',elename)
2226 open (ielep,file=elename,status='old',action='read')
2227 c print *,"Processor",myrank," opened file IELEP"
2228 call getenv_loc('SIDEPAR',sidename)
2229 open (isidep,file=sidename,status='old',action='read')
2230 c print *,"Processor",myrank," opened file ISIDEP"
2231 c print *,"Processor",myrank," opened parameter files"
2233 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2234 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2235 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2236 C Get parameter filenames and open the parameter files.
2237 call getenv_loc('BONDPAR',bondname)
2238 open (ibond,file=bondname,status='old')
2239 call getenv_loc('THETPAR',thetname)
2240 open (ithep,file=thetname,status='old')
2242 call getenv_loc('THETPARPDB',thetname_pdb)
2243 open (ithep_pdb,file=thetname_pdb,status='old')
2245 call getenv_loc('ROTPAR',rotname)
2246 open (irotam,file=rotname,status='old')
2248 call getenv_loc('ROTPARPDB',rotname_pdb)
2249 open (irotam_pdb,file=rotname_pdb,status='old')
2251 call getenv_loc('TORPAR',torname)
2252 open (itorp,file=torname,status='old')
2253 call getenv_loc('TORDPAR',tordname)
2254 open (itordp,file=tordname,status='old')
2255 call getenv_loc('SCCORPAR',sccorname)
2256 open (isccor,file=sccorname,status='old')
2257 call getenv_loc('FOURIER',fouriername)
2258 open (ifourier,file=fouriername,status='old')
2259 call getenv_loc('ELEPAR',elename)
2260 open (ielep,file=elename,status='old')
2261 call getenv_loc('SIDEPAR',sidename)
2262 open (isidep,file=sidename,status='old')
2264 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2266 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2267 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2268 C Get parameter filenames and open the parameter files.
2269 call getenv_loc('BONDPAR',bondname)
2270 open (ibond,file=bondname,status='old',action='read')
2271 call getenv_loc('THETPAR',thetname)
2272 open (ithep,file=thetname,status='old',action='read')
2274 call getenv_loc('THETPARPDB',thetname_pdb)
2275 print *,"thetname_pdb ",thetname_pdb
2276 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2277 print *,ithep_pdb," opened"
2279 call getenv_loc('ROTPAR',rotname)
2280 open (irotam,file=rotname,status='old',action='read')
2282 call getenv_loc('ROTPARPDB',rotname_pdb)
2283 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2285 call getenv_loc('TORPAR',torname)
2286 open (itorp,file=torname,status='old',action='read')
2287 call getenv_loc('TORDPAR',tordname)
2288 open (itordp,file=tordname,status='old',action='read')
2289 call getenv_loc('SCCORPAR',sccorname)
2290 open (isccor,file=sccorname,status='old',action='read')
2291 call getenv_loc('FOURIER',fouriername)
2292 open (ifourier,file=fouriername,status='old',action='read')
2293 call getenv_loc('ELEPAR',elename)
2294 open (ielep,file=elename,status='old',action='read')
2295 call getenv_loc('SIDEPAR',sidename)
2296 open (isidep,file=sidename,status='old',action='read')
2300 C 8/9/01 In the newest version SCp interaction constants are read from a file
2301 C Use -DOLDSCP to use hard-coded constants instead.
2303 call getenv_loc('SCPPAR',scpname)
2304 #if defined(WINIFL) || defined(WINPGI)
2305 open (iscpp,file=scpname,status='old',readonly,shared)
2306 #elif (defined CRAY) || (defined AIX)
2307 open (iscpp,file=scpname,status='old',action='read')
2309 open (iscpp,file=scpname,status='old')
2311 open (iscpp,file=scpname,status='old',action='read')
2314 call getenv_loc('PATTERN',patname)
2315 #if defined(WINIFL) || defined(WINPGI)
2316 open (icbase,file=patname,status='old',readonly,shared)
2317 #elif (defined CRAY) || (defined AIX)
2318 open (icbase,file=patname,status='old',action='read')
2320 open (icbase,file=patname,status='old')
2322 open (icbase,file=patname,status='old',action='read')
2325 C Open output file only for CG processes
2326 c print *,"Processor",myrank," fg_rank",fg_rank
2327 if (fg_rank.eq.0) then
2329 if (nodes.eq.1) then
2332 npos = dlog10(dfloat(nodes-1))+1
2334 if (npos.lt.3) npos=3
2335 write (liczba,'(i1)') npos
2336 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2338 write (liczba,form) me
2339 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2340 & liczba(:ilen(liczba))
2341 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2343 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2345 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2346 & liczba(:ilen(liczba))//'.mol2'
2347 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
2348 & liczba(:ilen(liczba))//'.x'
2349 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2350 & liczba(:ilen(liczba))//'.stat'
2352 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2353 & //liczba(:ilen(liczba))//'.stat')
2354 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2357 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2358 & liczba(:ilen(liczba))//'.const'
2363 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2364 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2365 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2366 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2367 cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
2368 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2370 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2372 rest2name=prefix(:ilen(prefix))//'.rst'
2374 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2377 #if defined(AIX) || defined(PGI)
2378 if (me.eq.king .or. .not. out1file)
2379 & open(iout,file=outname,status='unknown')
2382 if (fg_rank.gt.0) then
2383 write (liczba,'(i3.3)') myrank/nfgtasks
2384 write (ll,'(bz,i3.3)') fg_rank
2385 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2391 open(igeom,file=intname,status='unknown',position='append')
2392 open(ipdb,file=pdbname,status='unknown')
2393 open(imol2,file=mol2name,status='unknown')
2394 open(istat,file=statname,status='unknown',position='append')
2396 c1out open(iout,file=outname,status='unknown')
2399 if (me.eq.king .or. .not.out1file)
2400 & open(iout,file=outname,status='unknown')
2403 if (fg_rank.gt.0) then
2404 print "Processor",fg_rank," opening output file"
2405 write (liczba,'(i3.3)') myrank/nfgtasks
2406 write (ll,'(bz,i3.3)') fg_rank
2407 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2413 open(igeom,file=intname,status='unknown',access='append')
2414 open(ipdb,file=pdbname,status='unknown')
2415 open(imol2,file=mol2name,status='unknown')
2416 open(istat,file=statname,status='unknown',access='append')
2418 c1out open(iout,file=outname,status='unknown')
2421 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2422 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2423 csa csa_history=prefix(:lenpre)//'.CSA.history'
2424 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2425 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2426 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2427 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2428 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2429 csa csa_int=prefix(:lenpre)//'.int'
2430 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2431 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2432 csa csa_in=prefix(:lenpre)//'.CSA.in'
2433 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2436 write (iout,'(80(1h-))')
2437 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2438 write (iout,'(80(1h-))')
2439 write (iout,*) "Input file : ",
2440 & pref_orig(:ilen(pref_orig))//'.inp'
2441 write (iout,*) "Output file : ",
2442 & outname(:ilen(outname))
2444 write (iout,*) "Sidechain potential file : ",
2445 & sidename(:ilen(sidename))
2447 write (iout,*) "SCp potential file : ",
2448 & scpname(:ilen(scpname))
2450 write (iout,*) "Electrostatic potential file : ",
2451 & elename(:ilen(elename))
2452 write (iout,*) "Cumulant coefficient file : ",
2453 & fouriername(:ilen(fouriername))
2454 write (iout,*) "Torsional parameter file : ",
2455 & torname(:ilen(torname))
2456 write (iout,*) "Double torsional parameter file : ",
2457 & tordname(:ilen(tordname))
2458 write (iout,*) "SCCOR parameter file : ",
2459 & sccorname(:ilen(sccorname))
2460 write (iout,*) "Bond & inertia constant file : ",
2461 & bondname(:ilen(bondname))
2462 write (iout,*) "Bending parameter file : ",
2463 & thetname(:ilen(thetname))
2464 write (iout,*) "Rotamer parameter file : ",
2465 & rotname(:ilen(rotname))
2466 write (iout,*) "Threading database : ",
2467 & patname(:ilen(patname))
2469 &write (iout,*)" DIRTMP : ",
2471 write (iout,'(80(1h-))')
2475 c----------------------------------------------------------------------------
2476 subroutine card_concat(card)
2477 implicit real*8 (a-h,o-z)
2478 include 'DIMENSIONS'
2479 include 'COMMON.IOUNITS'
2481 character*80 karta,ucase
2483 read (inp,'(a)') karta
2486 do while (karta(80:80).eq.'&')
2487 card=card(:ilen(card)+1)//karta(:79)
2488 read (inp,'(a)') karta
2491 card=card(:ilen(card)+1)//karta
2494 c----------------------------------------------------------------------------------
2496 implicit real*8 (a-h,o-z)
2497 include 'DIMENSIONS'
2498 include 'COMMON.CHAIN'
2499 include 'COMMON.IOUNITS'
2501 include 'COMMON.CONTROL'
2502 open(irest2,file=rest2name,status='unknown')
2503 read(irest2,*) totT,EK,potE,totE,t_bath
2505 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2508 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2510 if(usampl.or.homol_nset.gt.1) then
2511 read (irest2,*) iset
2516 c---------------------------------------------------------------------------------
2517 subroutine read_fragments
2518 implicit real*8 (a-h,o-z)
2519 include 'DIMENSIONS'
2523 include 'COMMON.SETUP'
2524 include 'COMMON.CHAIN'
2525 include 'COMMON.IOUNITS'
2527 include 'COMMON.CONTROL'
2529 read(inp,*) nset,nfrag,npair,nfrag_back
2530 if(me.eq.king.or..not.out1file)
2531 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2532 & " nfrag_back",nfrag_back
2534 read(inp,*) mset(iset1)
2536 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2538 if(me.eq.king.or..not.out1file)
2539 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2540 & ifrag(2,i,iset1), qinfrag(i,iset1)
2543 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2545 if(me.eq.king.or..not.out1file)
2546 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2547 & ipair(2,i,iset1), qinpair(i,iset1)
2550 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2551 & wfrag_back(3,i,iset1),
2552 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2553 if(me.eq.king.or..not.out1file)
2554 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2555 & wfrag_back(2,i,iset1),
2556 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2557 & ifrag_back(2,i,iset1)
2562 c-------------------------------------------------------------------------------
2563 subroutine read_dist_constr
2564 implicit real*8 (a-h,o-z)
2565 include 'DIMENSIONS'
2569 include 'COMMON.SETUP'
2570 include 'COMMON.CONTROL'
2571 include 'COMMON.CHAIN'
2572 include 'COMMON.IOUNITS'
2573 include 'COMMON.SBRIDGE'
2574 integer ifrag_(2,100),ipair_(2,100)
2575 double precision wfrag_(100),wpair_(100)
2576 character*500 controlcard
2577 c write (iout,*) "Calling read_dist_constr"
2578 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2580 call card_concat(controlcard)
2581 call readi(controlcard,"NFRAG",nfrag_,0)
2582 call readi(controlcard,"NPAIR",npair_,0)
2583 call readi(controlcard,"NDIST",ndist_,0)
2584 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2586 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2587 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2588 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2589 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2590 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2591 c write (iout,*) "IFRAG"
2593 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2595 c write (iout,*) "IPAIR"
2597 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2599 if (.not.refstr .and. nfrag.gt.0) then
2601 & "ERROR: no reference structure to compute distance restraints"
2603 & "Restraints must be specified explicitly (NDIST=number)"
2606 if (nfrag.lt.2 .and. npair.gt.0) then
2607 write (iout,*) "ERROR: Less than 2 fragments specified",
2608 & " but distance restraints between pairs requested"
2613 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2614 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2615 & ifrag_(2,i)=nstart_sup+nsup-1
2616 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2618 if (wfrag_(i).gt.0.0d0) then
2619 do j=ifrag_(1,i),ifrag_(2,i)-1
2620 do k=j+1,ifrag_(2,i)
2621 c write (iout,*) "j",j," k",k
2623 if (constr_dist.eq.1) then
2628 forcon(nhpb)=wfrag_(i)
2629 else if (constr_dist.eq.2) then
2630 if (ddjk.le.dist_cut) then
2635 forcon(nhpb)=wfrag_(i)
2642 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2645 if (.not.out1file .or. me.eq.king)
2646 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2647 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2649 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2650 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2657 if (wpair_(i).gt.0.0d0) then
2665 do j=ifrag_(1,ii),ifrag_(2,ii)
2666 do k=ifrag_(1,jj),ifrag_(2,jj)
2670 forcon(nhpb)=wpair_(i)
2671 dhpb(nhpb)=dist(j,k)
2673 if (.not.out1file .or. me.eq.king)
2674 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2675 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2677 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2678 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2685 if (constr_dist.eq.11) then
2686 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2687 & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2688 fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2690 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2691 & ibecarb(i),forcon(nhpb+1)
2693 if (forcon(nhpb+1).gt.0.0d0) then
2695 if (ibecarb(i).gt.0) then
2696 ihpb(i)=ihpb(i)+nres
2697 jhpb(i)=jhpb(i)+nres
2699 if (dhpb(nhpb).eq.0.0d0)
2700 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2704 if (.not.out1file .or. me.eq.king) then
2707 if (constr_dist.eq.11) then
2708 write (iout,'(a,3i5,2f8.2,i2,2f10.1)') "+dist.constr11 ",
2709 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i),
2712 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2713 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2722 c-------------------------------------------------------------------------------
2724 subroutine read_constr_homology
2726 include 'DIMENSIONS'
2730 include 'COMMON.SETUP'
2731 include 'COMMON.CONTROL'
2732 include 'COMMON.CHAIN'
2733 include 'COMMON.IOUNITS'
2735 include 'COMMON.GEO'
2736 include 'COMMON.INTERACT'
2737 include 'COMMON.NAMES'
2739 c For new homol impl
2741 include 'COMMON.VAR'
2744 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2746 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2747 c & sigma_odl_temp(maxres,maxres,max_template)
2749 character*24 model_ki_dist, model_ki_angle
2750 character*500 controlcard
2751 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2752 logical lprn /.true./
2757 c FP - Nov. 2014 Temporary specifications for new vars
2759 double precision rescore_tmp,x12,y12,z12,rescore2_tmp
2760 double precision, dimension (max_template,maxres) :: rescore
2761 double precision, dimension (max_template,maxres) :: rescore2
2762 character*24 pdbfile,tpl_k_rescore
2763 c -----------------------------------------------------------------
2764 c Reading multiple PDB ref structures and calculation of retraints
2765 c not using pre-computed ones stored in files model_ki_{dist,angle}
2767 c -----------------------------------------------------------------
2770 c Alternative: reading from input
2771 call card_concat(controlcard)
2772 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2773 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2774 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2775 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2776 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2777 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2778 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2779 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2780 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2781 if(.not.read2sigma.and.start_from_model) then
2782 write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2783 start_from_model=.false.
2785 if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON'
2786 if(start_from_model .and. rest) then
2787 write(iout,*) 'START_FROM_MODELS is OFF'
2788 write(iout,*) 'remove restart keyword from input'
2790 if (homol_nset.gt.1)then
2791 call card_concat(controlcard)
2792 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2793 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2794 write(iout,*) "iset homology_weight "
2796 write(iout,*) i,waga_homology(i)
2799 iset=mod(kolor,homol_nset)+1
2802 waga_homology(1)=1.0
2805 cd write (iout,*) "nnt",nnt," nct",nct
2812 write(iout,*) 'nnt=',nnt,'nct=',nct
2815 do k=1,constr_homology
2828 do k=1,constr_homology
2830 read(inp,'(a)') pdbfile
2831 c Next stament causes error upon compilation (?)
2832 c if(me.eq.king.or. .not. out1file)
2833 c write (iout,'(2a)') 'PDB data will be read from file ',
2834 c & pdbfile(:ilen(pdbfile))
2835 write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2836 & pdbfile(:ilen(pdbfile))
2837 open(ipdbin,file=pdbfile,status='old',err=33)
2839 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2840 & pdbfile(:ilen(pdbfile))
2843 c print *,'Begin reading pdb data'
2845 c Files containing res sim or local scores (former containing sigmas)
2848 write(kic2,'(bz,i2.2)') k
2850 tpl_k_rescore="template"//kic2//".sco"
2853 if (read2sigma) then
2854 call readpdb_template(k)
2859 c Distance restraints
2862 C Copy the coordinates from reference coordinates (?)
2866 c write (iout,*) "c(",j,i,") =",c(j,i)
2870 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2872 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2873 open (ientin,file=tpl_k_rescore,status='old')
2874 if (nnt.gt.1) rescore(k,1)=0.0d0
2875 do irec=nnt,nct ! loop for reading res sim
2876 if (read2sigma) then
2877 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2880 idomain(k,i_tmp)=idomain_tmp
2881 rescore(k,i_tmp)=rescore_tmp
2882 rescore2(k,i_tmp)=rescore2_tmp
2883 write(iout,'(a7,i5,2f10.5,i5)') "rescore",
2884 & i_tmp,rescore2_tmp,rescore_tmp,
2888 read (ientin,*,end=1401) rescore_tmp
2890 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2891 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2892 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2897 if (waga_dist.ne.0.0d0) then
2905 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2906 c write (iout,*) k,i,j,distal,dist2_cut
2908 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2909 & .and. distal.le.dist2_cut ) then
2915 c write (iout,*) "k",k
2916 c write (iout,*) "i",i," j",j," constr_homology",
2921 if (read2sigma) then
2924 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2926 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2927 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2928 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2930 if (odl(k,ii).le.dist_cut) then
2931 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2934 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2935 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2937 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2938 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2942 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2945 l_homo(k,ii)=.false.
2952 c Theta, dihedral and SC retraints
2954 if (waga_angle.gt.0.0d0) then
2955 c open (ientin,file=tpl_k_sigma_dih,status='old')
2956 c do irec=1,maxres-3 ! loop for reading sigma_dih
2957 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2958 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2959 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2960 c & sigma_dih(k,i+nnt-1)
2965 if (idomain(k,i).eq.0) then
2969 dih(k,i)=phiref(i) ! right?
2970 c read (ientin,*) sigma_dih(k,i) ! original variant
2971 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2972 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2973 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2974 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2975 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2977 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2978 & rescore(k,i-2)+rescore(k,i-3))/4.0
2979 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2980 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2981 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2982 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2983 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2984 c Instead of res sim other local measure of b/b str reliability possible
2985 if (sigma_dih(k,i).ne.0)
2986 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2987 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2992 if (waga_theta.gt.0.0d0) then
2993 c open (ientin,file=tpl_k_sigma_theta,status='old')
2994 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2995 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2996 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2997 c & sigma_theta(k,i+nnt-1)
3002 do i = nnt+2,nct ! right? without parallel.
3003 c do i = i=1,nres ! alternative for bounds acc to readpdb?
3004 c do i=ithet_start,ithet_end ! with FG parallel.
3005 if (idomain(k,i).eq.0) then
3006 sigma_theta(k,i)=0.0
3009 thetatpl(k,i)=thetaref(i)
3010 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
3011 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
3012 c & "rescore(",k,i-1,") =",rescore(k,i-1),
3013 c & "rescore(",k,i-2,") =",rescore(k,i-2)
3014 c read (ientin,*) sigma_theta(k,i) ! 1st variant
3015 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
3016 & rescore(k,i-2))/3.0
3017 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
3018 if (sigma_theta(k,i).ne.0)
3019 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
3021 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
3022 c rescore(k,i-2) ! right expression ?
3023 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
3027 if (waga_d.gt.0.0d0) then
3028 c open (ientin,file=tpl_k_sigma_d,status='old')
3029 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3030 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3031 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3032 c & sigma_d(k,i+nnt-1)
3036 do i = nnt,nct ! right? without parallel.
3037 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
3038 c do i=loc_start,loc_end ! with FG parallel.
3039 if (itype(i).eq.10) cycle
3040 if (idomain(k,i).eq.0 ) then
3047 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3048 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3049 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3050 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3051 sigma_d(k,i)=rescore(k,i) ! right expression ?
3052 if (sigma_d(k,i).ne.0)
3053 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3055 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3056 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3057 c read (ientin,*) sigma_d(k,i) ! 1st variant
3058 if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
3063 c remove distance restraints not used in any model from the list
3064 c shift data in all arrays
3066 if (waga_dist.ne.0.0d0) then
3072 if (ii_in_use(ii).eq.0.and.liiflag) then
3076 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3077 & .not.liiflag.and.ii.eq.lim_odl) then
3078 if (ii.eq.lim_odl) then
3079 iishift=ii-iistart+1
3084 do ki=iistart,lim_odl-iishift
3085 ires_homo(ki)=ires_homo(ki+iishift)
3086 jres_homo(ki)=jres_homo(ki+iishift)
3087 ii_in_use(ki)=ii_in_use(ki+iishift)
3088 do k=1,constr_homology
3089 odl(k,ki)=odl(k,ki+iishift)
3090 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3091 l_homo(k,ki)=l_homo(k,ki+iishift)
3095 lim_odl=lim_odl-iishift
3100 if (constr_homology.gt.0) call homology_partition
3101 if (constr_homology.gt.0) call init_int_table
3102 cd write (iout,*) "homology_partition: lim_theta= ",lim_theta,
3103 cd & "lim_xx=",lim_xx
3104 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3105 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3109 if (.not.lprn) return
3110 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3111 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3112 write (iout,*) "Distance restraints from templates"
3114 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3115 & ii,ires_homo(ii),jres_homo(ii),
3116 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3117 & ki=1,constr_homology)
3119 write (iout,*) "Dihedral angle restraints from templates"
3121 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3122 & (rad2deg*dih(ki,i),
3123 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3125 write (iout,*) "Virtual-bond angle restraints from templates"
3127 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3128 & (rad2deg*thetatpl(ki,i),
3129 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3131 write (iout,*) "SC restraints from templates"
3133 write(iout,'(i5,100(4f8.2,4x))') i,
3134 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3135 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3138 c -----------------------------------------------------------------
3141 c----------------------------------------------------------------------
3144 subroutine flush(iu)
3149 subroutine flush(iu)
3155 c------------------------------------------------------------------------------
3156 subroutine copy_to_tmp(source)
3157 include "DIMENSIONS"
3158 include "COMMON.IOUNITS"
3159 character*(*) source
3160 character* 256 tmpfile
3164 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3165 inquire(file=tmpfile,exist=ex)
3167 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3168 & " to temporary directory..."
3169 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3170 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3174 c------------------------------------------------------------------------------
3175 subroutine move_from_tmp(source)
3176 include "DIMENSIONS"
3177 include "COMMON.IOUNITS"
3178 character*(*) source
3181 write (*,*) "Moving ",source(:ilen(source)),
3182 & " from temporary directory to working directory"
3183 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3184 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3187 c------------------------------------------------------------------------------
3188 subroutine random_init(seed)
3190 C Initialize random number generator
3192 implicit real*8 (a-h,o-z)
3193 include 'DIMENSIONS'
3199 logical OKRandom, prng_restart
3201 integer iseed_array(4)
3203 include 'COMMON.IOUNITS'
3204 include 'COMMON.TIME1'
3205 include 'COMMON.THREAD'
3206 include 'COMMON.SBRIDGE'
3207 include 'COMMON.CONTROL'
3208 include 'COMMON.MCM'
3209 include 'COMMON.MAP'
3210 include 'COMMON.HEADER'
3211 csa include 'COMMON.CSA'
3212 include 'COMMON.CHAIN'
3213 include 'COMMON.MUCA'
3215 include 'COMMON.FFIELD'
3216 include 'COMMON.SETUP'
3217 iseed=-dint(dabs(seed))
3218 if (iseed.eq.0) then
3219 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3220 & 'Random seed undefined. The program will stop.'
3221 write (*,'(/80(1h*)/20x,a/80(1h*))')
3222 & 'Random seed undefined. The program will stop.'
3224 call mpi_finalize(mpi_comm_world,ierr)
3226 stop 'Bad random seed.'
3229 if (fg_rank.eq.0) then
3233 if(me.eq.king .or. .not. out1file)
3234 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3235 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3236 OKRandom = prng_restart(me,iseedi8)
3239 tmp=65536.0d0**(4-i)
3240 iseed_array(i) = dint(seed/tmp)
3241 seed=seed-iseed_array(i)*tmp
3243 if(me.eq.king .or. .not. out1file)
3244 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3245 & (iseed_array(i),i=1,4)
3246 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3247 & (iseed_array(i),i=1,4)
3248 OKRandom = prng_restart(me,iseed_array)
3251 r1=ran_number(0.0D0,1.0D0)
3252 if(me.eq.king .or. .not. out1file)
3253 & write (iout,*) 'ran_num',r1
3254 if (r1.lt.0.0d0) OKRandom=.false.
3256 if (.not.OKRandom) then
3257 write (iout,*) 'PRNG IS NOT WORKING!!!'
3258 print *,'PRNG IS NOT WORKING!!!'
3261 call mpi_abort(mpi_comm_world,error_msg,ierr)
3264 write (iout,*) 'too many processors for parallel prng'
3265 write (*,*) 'too many processors for parallel prng'
3273 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)