2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
12 C Read force-field parameters except weights
14 C Read job setup parameters
16 C Read control parameters for energy minimzation if required
17 if (minim) call read_minim
18 C Read MCM control parameters if required
19 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
20 C Read MD control parameters if reqjuired
21 if (modecalc.eq.12) call read_MDpar
22 C Read MREMD control parameters if required
23 if (modecalc.eq.14) then
27 C Read MUCA control parameters if required
28 if (lmuca) call read_muca
29 C Read CSA control parameters if required (from fort.40 if exists
30 C otherwise from general input file)
31 csa if (modecalc.eq.8) then
32 csa inquire (file="fort.40",exist=file_exist)
33 csa if (.not.file_exist) call csaread
35 cfmc if (modecalc.eq.10) call mcmfread
36 C Read molecule information, molecule geometry, energy-term weights, and
37 C restraints if requested
39 C Print restraint information
41 if (.not. out1file .or. me.eq.king) then
44 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
45 & " distance constraints have been imposed"
47 write (iout,'(3i6,i2,3f10.5)') i-nss,ihpb(i),jhpb(i),
48 & ibecarb(i),dhpb(i),dhpb1(i),forcon(i)
53 c print *,"Processor",myrank," leaves READRTNS"
56 C-------------------------------------------------------------------------------
57 subroutine read_control
61 implicit real*8 (a-h,o-z)
65 logical OKRandom, prng_restart
68 include 'COMMON.IOUNITS'
69 include 'COMMON.TIME1'
70 include 'COMMON.THREAD'
71 include 'COMMON.SBRIDGE'
72 include 'COMMON.CONTROL'
75 include 'COMMON.HEADER'
76 csa include 'COMMON.CSA'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.SETUP'
82 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
83 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
85 character*320 controlcard
90 read (INP,'(a)') titel
91 call card_concat(controlcard)
92 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
93 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
94 call reada(controlcard,'SEED',seed,0.0D0)
95 call random_init(seed)
96 C Set up the time limit (caution! The time must be input in minutes!)
97 read_cart=index(controlcard,'READ_CART').gt.0
98 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
99 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
100 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
101 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
102 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
103 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
104 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
105 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
106 call reada(controlcard,'DRMS',drms,0.1D0)
107 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
108 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
109 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
110 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
111 write (iout,'(a,f10.1)')'DRMS = ',drms
112 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
113 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
115 call readi(controlcard,'NZ_START',nz_start,0)
116 call readi(controlcard,'NZ_END',nz_end,0)
117 call readi(controlcard,'IZ_SC',iz_sc,0)
119 safety = 60.0d0*safety
122 call reada(controlcard,"T_BATH",t_bath,300.0d0)
123 minim=(index(controlcard,'MINIMIZE').gt.0)
124 dccart=(index(controlcard,'CART').gt.0)
125 overlapsc=(index(controlcard,'OVERLAP').gt.0)
126 overlapsc=.not.overlapsc
127 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
128 searchsc=.not.searchsc
129 sideadd=(index(controlcard,'SIDEADD').gt.0)
130 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
131 outpdb=(index(controlcard,'PDBOUT').gt.0)
132 outmol2=(index(controlcard,'MOL2OUT').gt.0)
133 pdbref=(index(controlcard,'PDBREF').gt.0)
134 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
135 indpdb=index(controlcard,'PDBSTART')
136 extconf=(index(controlcard,'EXTCONF').gt.0)
137 call readi(controlcard,'IPRINT',iprint,0)
138 call readi(controlcard,'MAXGEN',maxgen,10000)
139 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
140 call readi(controlcard,"KDIAG",kdiag,0)
141 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
142 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
143 & write (iout,*) "RESCALE_MODE",rescale_mode
144 split_ene=index(controlcard,'SPLIT_ENE').gt.0
145 if (index(controlcard,'REGULAR').gt.0.0D0) then
146 call reada(controlcard,'WEIDIS',weidis,0.1D0)
150 if (index(controlcard,'CHECKGRAD').gt.0) then
152 if (index(controlcard,'CART').gt.0) then
154 elseif (index(controlcard,'CARINT').gt.0) then
159 elseif (index(controlcard,'THREAD').gt.0) then
161 call readi(controlcard,'THREAD',nthread,0)
162 if (nthread.gt.0) then
163 call reada(controlcard,'WEIDIS',weidis,0.1D0)
166 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
167 stop 'Error termination in Read_Control.'
169 else if (index(controlcard,'MCMA').gt.0) then
171 else if (index(controlcard,'MCEE').gt.0) then
173 else if (index(controlcard,'MULTCONF').gt.0) then
175 else if (index(controlcard,'MAP').gt.0) then
177 call readi(controlcard,'MAP',nmap,0)
178 else if (index(controlcard,'CSA').gt.0) then
179 write(*,*) "CSA not supported in this version"
182 crc else if (index(controlcard,'ZSCORE').gt.0) then
184 crc ZSCORE is rm from UNRES, modecalc=9 is available
187 cfcm else if (index(controlcard,'MCMF').gt.0) then
189 else if (index(controlcard,'SOFTREG').gt.0) then
191 else if (index(controlcard,'CHECK_BOND').gt.0) then
193 else if (index(controlcard,'TEST').gt.0) then
195 else if (index(controlcard,'MD').gt.0) then
197 else if (index(controlcard,'RE ').gt.0) then
201 lmuca=index(controlcard,'MUCA').gt.0
202 call readi(controlcard,'MUCADYN',mucadyn,0)
203 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
204 if (lmuca .and. (me.eq.king .or. .not.out1file ))
206 write (iout,*) 'MUCADYN=',mucadyn
207 write (iout,*) 'MUCASMOOTH=',muca_smooth
210 iscode=index(controlcard,'ONE_LETTER')
211 indphi=index(controlcard,'PHI')
212 indback=index(controlcard,'BACK')
213 iranconf=index(controlcard,'RAND_CONF')
214 i2ndstr=index(controlcard,'USE_SEC_PRED')
215 gradout=index(controlcard,'GRADOUT').gt.0
216 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
218 if(me.eq.king.or..not.out1file)
219 & write (iout,'(2a)') diagmeth(kdiag),
220 & ' routine used to diagonalize matrices.'
223 c--------------------------------------------------------------------------
224 subroutine read_REMDpar
228 implicit real*8 (a-h,o-z)
230 include 'COMMON.IOUNITS'
231 include 'COMMON.TIME1'
234 include 'COMMON.LANGEVIN'
236 include 'COMMON.LANGEVIN.lang0'
238 include 'COMMON.INTERACT'
239 include 'COMMON.NAMES'
241 include 'COMMON.REMD'
242 include 'COMMON.CONTROL'
243 include 'COMMON.SETUP'
245 character*320 controlcard
246 character*3200 controlcard1
247 integer iremd_m_total
249 if(me.eq.king.or..not.out1file)
250 & write (iout,*) "REMD setup"
252 call card_concat(controlcard)
253 call readi(controlcard,"NREP",nrep,3)
254 call readi(controlcard,"NSTEX",nstex,1000)
255 call reada(controlcard,"RETMIN",retmin,10.0d0)
256 call reada(controlcard,"RETMAX",retmax,1000.0d0)
257 mremdsync=(index(controlcard,'SYNC').gt.0)
258 call readi(controlcard,"NSYN",i_sync_step,100)
259 restart1file=(index(controlcard,'REST1FILE').gt.0)
260 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
261 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
262 if(max_cache_traj_use.gt.max_cache_traj)
263 & max_cache_traj_use=max_cache_traj
264 if(me.eq.king.or..not.out1file) then
265 cd if (traj1file) then
266 crc caching is in testing - NTWX is not ignored
267 cd write (iout,*) "NTWX value is ignored"
268 cd write (iout,*) " trajectory is stored to one file by master"
269 cd write (iout,*) " before exchange at NSTEX intervals"
271 write (iout,*) "NREP= ",nrep
272 write (iout,*) "NSTEX= ",nstex
273 write (iout,*) "SYNC= ",mremdsync
274 write (iout,*) "NSYN= ",i_sync_step
275 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
278 t_exchange_only=(index(controlcard,'TONLY').gt.0)
279 call readi(controlcard,"HREMD",hremd,0)
280 if((me.eq.king.or..not.out1file).and.hremd.gt.0) then
281 write (iout,*) "Hamiltonian REMD with ",hremd," sets of weights"
283 if(usampl.and.hremd.gt.0) then
285 & "========== ERROR: USAMPL and HREMD cannot be used together"
287 call MPI_Abort(MPI_COMM_WORLD,IERROR,ERRCODE)
294 if (index(controlcard,'TLIST').gt.0) then
296 call card_concat(controlcard1)
297 read(controlcard1,*) (remd_t(i),i=1,nrep)
298 if(me.eq.king.or..not.out1file)
299 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
302 if (index(controlcard,'MLIST').gt.0) then
304 call card_concat(controlcard1)
305 read(controlcard1,*) (remd_m(i),i=1,nrep)
306 if(me.eq.king.or..not.out1file) then
307 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
310 iremd_m_total=iremd_m_total+remd_m(i)
313 write (iout,*) 'Total number of replicas ',
314 & iremd_m_total*hremd
316 write (iout,*) 'Total number of replicas ',iremd_m_total
320 if(me.eq.king.or..not.out1file)
321 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
324 c--------------------------------------------------------------------------
325 subroutine read_MDpar
329 implicit real*8 (a-h,o-z)
331 include 'COMMON.IOUNITS'
332 include 'COMMON.TIME1'
335 include 'COMMON.LANGEVIN'
337 include 'COMMON.LANGEVIN.lang0'
339 include 'COMMON.INTERACT'
340 include 'COMMON.NAMES'
342 include 'COMMON.SETUP'
343 include 'COMMON.CONTROL'
344 include 'COMMON.SPLITELE'
346 character*320 controlcard
348 call card_concat(controlcard)
349 call readi(controlcard,"NSTEP",n_timestep,1000000)
350 call readi(controlcard,"NTWE",ntwe,100)
351 call readi(controlcard,"NTWX",ntwx,1000)
352 call reada(controlcard,"DT",d_time,1.0d-1)
353 call reada(controlcard,"DVMAX",dvmax,2.0d1)
354 call reada(controlcard,"DAMAX",damax,1.0d1)
355 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
356 call readi(controlcard,"LANG",lang,0)
357 RESPA = index(controlcard,"RESPA") .gt. 0
358 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
359 ntime_split0=ntime_split
360 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
361 ntime_split0=ntime_split
362 call reada(controlcard,"R_CUT",r_cut,2.0d0)
363 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
364 rest = index(controlcard,"REST").gt.0
365 tbf = index(controlcard,"TBF").gt.0
366 call readi(controlcard,"HMC",hmc,0)
367 tnp = index(controlcard,"NOSEPOINCARE99").gt.0
368 tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0
369 tnh = index(controlcard,"NOSEHOOVER96").gt.0
370 if (RESPA.and.tnh)then
371 xiresp = index(controlcard,"XIRESP").gt.0
373 call reada(controlcard,"Q_NP",Q_np,0.1d0)
374 usampl = index(controlcard,"USAMPL").gt.0
376 mdpdb = index(controlcard,"MDPDB").gt.0
377 call reada(controlcard,"T_BATH",t_bath,300.0d0)
378 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
379 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
380 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
381 if (count_reset_moment.eq.0) count_reset_moment=1000000000
382 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
383 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
384 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
385 if (count_reset_vel.eq.0) count_reset_vel=1000000000
386 large = index(controlcard,"LARGE").gt.0
387 print_compon = index(controlcard,"PRINT_COMPON").gt.0
388 rattle = index(controlcard,"RATTLE").gt.0
389 c if performing umbrella sampling, fragments constrained are read from the fragment file
395 if(me.eq.king.or..not.out1file) then
397 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
399 write (iout,'(a)') "The units are:"
400 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
401 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
402 & " acceleration: angstrom/(48.9 fs)**2"
403 write (iout,'(a)') "energy: kcal/mol, temperature: K"
405 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
406 write (iout,'(a60,f10.5,a)')
407 & "Initial time step of numerical integration:",d_time,
409 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
411 write (iout,'(2a,i4,a)')
412 & "A-MTS algorithm used; initial time step for fast-varying",
413 & " short-range forces split into",ntime_split," steps."
414 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
415 & r_cut," lambda",rlamb
417 write (iout,'(2a,f10.5)')
418 & "Maximum acceleration threshold to reduce the time step",
419 & "/increase split number:",damax
420 write (iout,'(2a,f10.5)')
421 & "Maximum predicted energy drift to reduce the timestep",
422 & "/increase split number:",edriftmax
423 write (iout,'(a60,f10.5)')
424 & "Maximum velocity threshold to reduce velocities:",dvmax
425 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
426 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
427 if (rattle) write (iout,'(a60)')
428 & "Rattle algorithm used to constrain the virtual bonds"
432 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
433 call reada(controlcard,"RWAT",rwat,1.4d0)
434 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
435 surfarea=index(controlcard,"SURFAREA").gt.0
436 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
437 if(me.eq.king.or..not.out1file)then
438 write (iout,'(/a,$)') "Langevin dynamics calculation"
441 & " with direct integration of Langevin equations"
442 else if (lang.eq.2) then
443 write (iout,'(a/)') " with TINKER stochasic MD integrator"
444 else if (lang.eq.3) then
445 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
446 else if (lang.eq.4) then
447 write (iout,'(a/)') " in overdamped mode"
449 write (iout,'(//a,i5)')
450 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
453 write (iout,'(a60,f10.5)') "Temperature:",t_bath
454 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
455 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
456 write (iout,'(a60,f10.5)')
457 & "Scaling factor of the friction forces:",scal_fric
458 if (surfarea) write (iout,'(2a,i10,a)')
459 & "Friction coefficients will be scaled by solvent-accessible",
460 & " surface area every",reset_fricmat," steps."
462 c Calculate friction coefficients and bounds of stochastic forces
463 eta=6*pi*cPoise*etawat
464 if(me.eq.king.or..not.out1file)
465 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
467 gamp=scal_fric*(pstok+rwat)*eta
468 stdfp=dsqrt(2*Rb*t_bath/d_time)
470 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
471 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
473 if(me.eq.king.or..not.out1file)then
474 write (iout,'(/2a/)')
475 & "Radii of site types and friction coefficients and std's of",
476 & " stochastic forces of fully exposed sites"
477 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
479 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
480 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
484 if(me.eq.king.or..not.out1file)then
485 write (iout,'(a)') "Berendsen bath calculation"
486 write (iout,'(a60,f10.5)') "Temperature:",t_bath
487 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
489 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
490 & count_reset_moment," steps"
492 & write (iout,'(a,i10,a)')
493 & "Velocities will be reset at random every",count_reset_vel,
496 else if (tnp .or. tnp1 .or. tnh) then
497 if (tnp .or. tnp1) then
498 write (iout,'(a)') "Nose-Poincare bath calculation"
499 if (tnp) write (iout,'(a)')
500 & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird"
501 if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose"
503 write (iout,'(a)') "Nose-Hoover bath calculation"
504 write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al."
514 WDTI(i) = 1.0*d_time/nresn
521 write (iout,'(a)') "NVT-XI-RESPA algorithm"
523 write (iout,'(a)') "NVT-XO-RESPA algorithm"
526 WDTIi(i) = 1.0*d_time/nresn/ntime_split
534 write (iout,'(a60,f10.5)') "Temperature:",t_bath
535 write (iout,'(a60,f10.5)') "Q =",Q_np
537 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
538 & count_reset_moment," steps"
540 & write (iout,'(a,i10,a)')
541 & "Velocities will be reset at random every",count_reset_vel,
544 else if (hmc.gt.0) then
545 write (iout,'(a)') "Hybrid Monte Carlo calculation"
546 write (iout,'(a60,f10.5)') "Temperature:",t_bath
547 write (iout,'(a60,i10)')
548 & "Number of MD steps between Metropolis tests:",hmc
551 if(me.eq.king.or..not.out1file)
552 & write (iout,'(a31)') "Microcanonical mode calculation"
554 if(me.eq.king.or..not.out1file)then
555 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
557 write(iout,*) "MD running with constraints."
558 write(iout,*) "Equilibration time ", eq_time, " mtus."
559 write(iout,*) "Constraining ", nfrag," fragments."
560 write(iout,*) "Length of each fragment, weight and q0:"
562 write (iout,*) "Set of restraints #",iset
564 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
565 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
567 write(iout,*) "constraints between ", npair, "fragments."
568 write(iout,*) "constraint pairs, weights and q0:"
570 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
571 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
573 write(iout,*) "angle constraints within ", nfrag_back,
574 & "backbone fragments."
575 write(iout,*) "fragment, weights:"
577 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
578 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
579 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
582 iset=mod(kolor,nset)+1
585 if(me.eq.king.or..not.out1file)
586 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
589 c------------------------------------------------------------------------------
592 C Read molecular data.
594 implicit real*8 (a-h,o-z)
600 include 'COMMON.IOUNITS'
603 include 'COMMON.INTERACT'
604 include 'COMMON.LOCAL'
605 include 'COMMON.NAMES'
606 include 'COMMON.CHAIN'
607 include 'COMMON.FFIELD'
608 include 'COMMON.SBRIDGE'
609 include 'COMMON.HEADER'
610 include 'COMMON.CONTROL'
611 include 'COMMON.DBASE'
612 include 'COMMON.THREAD'
613 include 'COMMON.CONTACTS'
614 include 'COMMON.TORCNSTR'
615 include 'COMMON.TIME1'
616 include 'COMMON.BOUNDS'
618 include 'COMMON.REMD'
619 include 'COMMON.SETUP'
620 character*4 sequence(maxres)
622 double precision x(maxvar)
623 character*256 pdbfile
624 character*320 weightcard
625 character*80 weightcard_t,ucase
626 dimension itype_pdb(maxres)
627 common /pizda/ itype_pdb
628 logical seq_comp,fail
629 double precision energia(0:n_ene)
635 C Read weights of the subsequent energy terms.
648 if(me.eq.king.or..not.out1file) then
649 write (iout,*) 'Reading ',hremd,' sets of weights for HREMD'
650 write (iout,*) 'Current weights for processor ',
651 & me,' set ',i2set(me)
655 call card_concat(weightcard)
656 call reada(weightcard,'WLONG',wlong,1.0D0)
657 call reada(weightcard,'WSC',wsc,wlong)
658 call reada(weightcard,'WSCP',wscp,wlong)
659 call reada(weightcard,'WELEC',welec,1.0D0)
660 call reada(weightcard,'WVDWPP',wvdwpp,welec)
661 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
662 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
663 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
664 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
665 call reada(weightcard,'WTURN3',wturn3,1.0D0)
666 call reada(weightcard,'WTURN4',wturn4,1.0D0)
667 call reada(weightcard,'WTURN6',wturn6,1.0D0)
668 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
669 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
670 call reada(weightcard,'WBOND',wbond,1.0D0)
671 call reada(weightcard,'WTOR',wtor,1.0D0)
672 call reada(weightcard,'WTORD',wtor_d,1.0D0)
673 call reada(weightcard,'WANG',wang,1.0D0)
674 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
675 call reada(weightcard,'SCAL14',scal14,0.4D0)
676 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
677 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
678 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
679 call reada(weightcard,'TEMP0',temp0,300.0d0)
680 if (index(weightcard,'SOFT').gt.0) ipot=6
681 C 12/1/95 Added weight for the multi-body term WCORR
682 call reada(weightcard,'WCORRH',wcorr,1.0D0)
683 if (wcorr4.gt.0.0d0) wcorr=wcorr4
691 hweights(i,7)=wel_loc
694 hweights(i,10)=wturn6
696 hweights(i,12)=wscloc
698 hweights(i,14)=wtor_d
699 hweights(i,15)=wstrain
700 hweights(i,16)=wvdwpp
702 hweights(i,18)=scal14
703 hweights(i,21)=wsccor
708 weights(i)=hweights(i2set(me),i)
732 call card_concat(weightcard)
733 call reada(weightcard,'WLONG',wlong,1.0D0)
734 call reada(weightcard,'WSC',wsc,wlong)
735 call reada(weightcard,'WSCP',wscp,wlong)
736 call reada(weightcard,'WELEC',welec,1.0D0)
737 call reada(weightcard,'WVDWPP',wvdwpp,welec)
738 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
739 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
740 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
741 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
742 call reada(weightcard,'WTURN3',wturn3,1.0D0)
743 call reada(weightcard,'WTURN4',wturn4,1.0D0)
744 call reada(weightcard,'WTURN6',wturn6,1.0D0)
745 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
746 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
747 call reada(weightcard,'WBOND',wbond,1.0D0)
748 call reada(weightcard,'WTOR',wtor,1.0D0)
749 call reada(weightcard,'WTORD',wtor_d,1.0D0)
750 call reada(weightcard,'WANG',wang,1.0D0)
751 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
752 call reada(weightcard,'SCAL14',scal14,0.4D0)
753 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
754 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
755 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
756 call reada(weightcard,'TEMP0',temp0,300.0d0)
757 if (index(weightcard,'SOFT').gt.0) ipot=6
758 C 12/1/95 Added weight for the multi-body term WCORR
759 call reada(weightcard,'WCORRH',wcorr,1.0D0)
760 if (wcorr4.gt.0.0d0) wcorr=wcorr4
782 if(me.eq.king.or..not.out1file)
783 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
784 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
786 10 format (/'Energy-term weights (unscaled):'//
787 & 'WSCC= ',f10.6,' (SC-SC)'/
788 & 'WSCP= ',f10.6,' (SC-p)'/
789 & 'WELEC= ',f10.6,' (p-p electr)'/
790 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
791 & 'WBOND= ',f10.6,' (stretching)'/
792 & 'WANG= ',f10.6,' (bending)'/
793 & 'WSCLOC= ',f10.6,' (SC local)'/
794 & 'WTOR= ',f10.6,' (torsional)'/
795 & 'WTORD= ',f10.6,' (double torsional)'/
796 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
797 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
798 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
799 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
800 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
801 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
802 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
803 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
804 & 'WTURN6= ',f10.6,' (turns, 6th order)')
805 if(me.eq.king.or..not.out1file)then
806 if (wcorr4.gt.0.0d0) then
807 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
808 & 'between contact pairs of peptide groups'
809 write (iout,'(2(a,f5.3/))')
810 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
811 & 'Range of quenching the correlation terms:',2*delt_corr
812 else if (wcorr.gt.0.0d0) then
813 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
814 & 'between contact pairs of peptide groups'
816 write (iout,'(a,f8.3)')
817 & 'Scaling factor of 1,4 SC-p interactions:',scal14
818 write (iout,'(a,f8.3)')
819 & 'General scaling factor of SC-p interactions:',scalscp
821 r0_corr=cutoff_corr-delt_corr
823 aad(i,1)=scalscp*aad(i,1)
824 aad(i,2)=scalscp*aad(i,2)
825 bad(i,1)=scalscp*bad(i,1)
826 bad(i,2)=scalscp*bad(i,2)
828 call rescale_weights(t_bath)
829 if(me.eq.king.or..not.out1file)
830 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
831 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
833 22 format (/'Energy-term weights (scaled):'//
834 & 'WSCC= ',f10.6,' (SC-SC)'/
835 & 'WSCP= ',f10.6,' (SC-p)'/
836 & 'WELEC= ',f10.6,' (p-p electr)'/
837 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
838 & 'WBOND= ',f10.6,' (stretching)'/
839 & 'WANG= ',f10.6,' (bending)'/
840 & 'WSCLOC= ',f10.6,' (SC local)'/
841 & 'WTOR= ',f10.6,' (torsional)'/
842 & 'WTORD= ',f10.6,' (double torsional)'/
843 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
844 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
845 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
846 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
847 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
848 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
849 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
850 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
851 & 'WTURN6= ',f10.6,' (turns, 6th order)')
852 if(me.eq.king.or..not.out1file)
853 & write (iout,*) "Reference temperature for weights calculation:",
855 call reada(weightcard,"D0CM",d0cm,3.78d0)
856 call reada(weightcard,"AKCM",akcm,15.1d0)
857 call reada(weightcard,"AKTH",akth,11.0d0)
858 call reada(weightcard,"AKCT",akct,12.0d0)
859 call reada(weightcard,"V1SS",v1ss,-1.08d0)
860 call reada(weightcard,"V2SS",v2ss,7.61d0)
861 call reada(weightcard,"V3SS",v3ss,13.7d0)
862 call reada(weightcard,"EBR",ebr,-5.50D0)
863 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
865 dyn_ss_mask(i)=.false.
869 dyn_ssbond_ij(i,j)=1.0d300
872 call reada(weightcard,"HT",Ht,0.0D0)
874 ss_depth=ebr/wsc-0.25*eps(1,1)
875 Ht=Ht/wsc-0.25*eps(1,1)
876 akcm=akcm*wstrain/wsc
877 akth=akth*wstrain/wsc
878 akct=akct*wstrain/wsc
879 v1ss=v1ss*wstrain/wsc
880 v2ss=v2ss*wstrain/wsc
881 v3ss=v3ss*wstrain/wsc
883 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
886 if(me.eq.king.or..not.out1file) then
887 write (iout,*) "Parameters of the SS-bond potential:"
888 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
890 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
891 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
892 write (iout,*)" HT",Ht
893 print *,'indpdb=',indpdb,' pdbref=',pdbref
895 if (indpdb.gt.0 .or. pdbref) then
896 read(inp,'(a)') pdbfile
897 if(me.eq.king.or..not.out1file)
898 & write (iout,'(2a)') 'PDB data will be read from file ',
899 & pdbfile(:ilen(pdbfile))
900 open(ipdbin,file=pdbfile,status='old',err=33)
902 33 write (iout,'(a)') 'Error opening PDB file.'
905 c print *,'Begin reading pdb data'
907 c print *,'Finished reading pdb data'
908 if(me.eq.king.or..not.out1file)
909 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
910 & ' nstart_sup=',nstart_sup
912 itype_pdb(i)=itype(i)
916 nct=nstart_sup+nsup-1
917 call contact(.false.,ncont_ref,icont_ref,co)
920 if(me.eq.king.or..not.out1file)
921 & write(iout,*)'Adding sidechains'
928 do while (fail.and.nsi.le.maxsi)
929 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
932 if(fail) write(iout,*)'Adding sidechain failed for res ',
933 & i,' after ',nsi,' trials'
938 if (indpdb.eq.0) then
939 C Read sequence if not taken from the pdb file.
941 c print *,'nres=',nres
942 if (iscode.gt.0) then
943 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
945 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
947 C Convert sequence to numeric code
949 itype(i)=rescode(i,sequence(i),iscode)
951 C Assign initial virtual bond lengths
957 vbld(i+nres)=dsc(itype(i))
958 vbld_inv(i+nres)=dsc_inv(itype(i))
959 c write (iout,*) "i",i," itype",itype(i),
960 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
964 c print '(20i4)',(itype(i),i=1,nres)
967 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
969 if (itype(i).eq.21) then
973 else if (itype(i+1).ne.20) then
975 else if (itype(i).ne.20) then
982 if(me.eq.king.or..not.out1file)then
983 write (iout,*) "ITEL"
985 write (iout,*) i,itype(i),itel(i)
987 print *,'Call Read_Bridge.'
990 C 8/13/98 Set limits to generating the dihedral angles
995 read (inp,*) ndih_constr
996 if (ndih_constr.gt.0) then
998 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
999 if(me.eq.king.or..not.out1file)then
1001 & 'There are',ndih_constr,' constraints on phi angles.'
1003 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1007 phi0(i)=deg2rad*phi0(i)
1008 drange(i)=deg2rad*drange(i)
1010 if(me.eq.king.or..not.out1file)
1011 & write (iout,*) 'FTORS',ftors
1014 phibound(1,ii) = phi0(i)-drange(i)
1015 phibound(2,ii) = phi0(i)+drange(i)
1020 if (me.eq.king) then
1022 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1024 write (iout,'(a3,i5,2f10.1)')
1025 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1031 cd print *,'NNT=',NNT,' NCT=',NCT
1032 if (itype(1).eq.21) nnt=2
1033 if (itype(nres).eq.21) nct=nct-1
1035 if(me.eq.king.or..not.out1file)
1036 & write (iout,'(a,i3)') 'nsup=',nsup
1038 if (nsup.le.(nct-nnt+1)) then
1039 do i=0,nct-nnt+1-nsup
1040 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1046 & 'Error - sequences to be superposed do not match.'
1049 do i=0,nsup-(nct-nnt+1)
1050 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1052 nstart_sup=nstart_sup+i
1058 & 'Error - sequences to be superposed do not match.'
1061 if (nsup.eq.0) nsup=nct-nnt
1062 if (nstart_sup.eq.0) nstart_sup=nnt
1063 if (nstart_seq.eq.0) nstart_seq=nnt
1064 if(me.eq.king.or..not.out1file)
1065 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1066 & ' nstart_seq=',nstart_seq
1068 c--- Zscore rms -------
1069 if (nz_start.eq.0) nz_start=nnt
1070 if (nz_end.eq.0 .and. nsup.gt.0) then
1072 else if (nz_end.eq.0) then
1075 if(me.eq.king.or..not.out1file)then
1076 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1077 write (iout,*) 'IZ_SC=',iz_sc
1079 c----------------------
1082 if (.not.pdbref) then
1083 call read_angles(inp,*38)
1085 38 write (iout,'(a)') 'Error reading reference structure.'
1087 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1088 stop 'Error reading reference structure'
1092 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1101 call contact(.true.,ncont_ref,icont_ref,co)
1103 if(me.eq.king.or..not.out1file)
1104 & write (iout,*) 'Contact order:',co
1106 if(me.eq.king.or..not.out1file)
1107 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1110 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1112 if(me.eq.king.or..not.out1file)
1113 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1114 & icont_ref(1,i),' ',
1115 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1119 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1120 if (constr_dist.gt.0) then
1121 call read_dist_constr
1123 if (nhpb.gt.0) call hpb_partition
1124 c write (iout,*) "After read_dist_constr nhpb",nhpb
1126 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1127 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1128 & modecalc.ne.10) then
1129 C If input structure hasn't been supplied from the PDB file read or generate
1131 if (iranconf.eq.0 .and. .not. extconf) then
1132 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1133 & write (iout,'(a)') 'Initial geometry will be read in.'
1135 read(inp,'(8f10.5)',end=36,err=36)
1136 & ((c(l,k),l=1,3),k=1,nres),
1137 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1138 call int_from_cart1(.false.)
1141 dc(j,i)=c(j,i+1)-c(j,i)
1142 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1146 if (itype(i).ne.10) then
1148 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1149 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1155 call read_angles(inp,*36)
1158 36 write (iout,'(a)') 'Error reading angle file.'
1160 call mpi_finalize( MPI_COMM_WORLD,IERR )
1162 stop 'Error reading angle file.'
1164 else if (extconf) then
1165 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1166 & write (iout,'(a)') 'Extended chain initial geometry.'
1168 theta(i)=90d0*deg2rad
1171 phi(i)=180d0*deg2rad
1174 alph(i)=110d0*deg2rad
1177 omeg(i)=-120d0*deg2rad
1180 if(me.eq.king.or..not.out1file)
1181 & write (iout,'(a)') 'Random-generated initial geometry.'
1185 if (me.eq.king .or. fg_rank.eq.0 .and. (
1186 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1190 call gen_rand_conf(itmp,*30)
1192 30 write (iout,*) 'Failed to generate random conformation',
1193 & ', itrial=',itrial
1194 write (*,*) 'Processor:',me,
1195 & ' Failed to generate random conformation',
1205 write (iout,'(a,i3,a)') 'Processor:',me,
1206 & ' error in generating random conformation.'
1207 write (*,'(a,i3,a)') 'Processor:',me,
1208 & ' error in generating random conformation.'
1211 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1218 elseif (modecalc.eq.4) then
1219 read (inp,'(a)') intinname
1220 open (intin,file=intinname,status='old',err=333)
1221 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1222 & write (iout,'(a)') 'intinname',intinname
1223 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1225 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1227 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1229 stop 'Error opening angle file.'
1233 C Generate distance constraints, if the PDB structure is to be regularized.
1234 if (nthread.gt.0) then
1235 call read_threadbase
1238 if (me.eq.king .or. .not. out1file)
1240 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1241 write (iout,'(/a,i3,a)')
1242 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1243 write (iout,'(20i4)') (iss(i),i=1,ns)
1245 write(iout,*)"Running with dynamic disulfide-bond formation"
1249 forcon(i-nss)=forcon(i)
1256 dyn_ss_mask(iss(i))=.true.
1259 write (iout,'(/a/)') 'Pre-formed links are:'
1265 if (me.eq.king.or..not.out1file)
1266 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1267 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1273 if (i2ndstr.gt.0) call secstrp2dihc
1274 c call geom_to_var(nvar,x)
1275 c call etotal(energia(0))
1276 c call enerprint(energia(0))
1277 c call briefout(0,etot)
1279 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1280 cd write (iout,'(a)') 'Variable list:'
1281 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1283 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1284 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1285 & 'Processor',myrank,': end reading molecular data.'
1289 c--------------------------------------------------------------------------
1290 logical function seq_comp(itypea,itypeb,length)
1292 integer length,itypea(length),itypeb(length)
1295 if (itypea(i).ne.itypeb(i)) then
1303 c-----------------------------------------------------------------------------
1304 subroutine read_bridge
1305 C Read information about disulfide bridges.
1306 implicit real*8 (a-h,o-z)
1307 include 'DIMENSIONS'
1311 include 'COMMON.IOUNITS'
1312 include 'COMMON.GEO'
1313 include 'COMMON.VAR'
1314 include 'COMMON.INTERACT'
1315 include 'COMMON.LOCAL'
1316 include 'COMMON.NAMES'
1317 include 'COMMON.CHAIN'
1318 include 'COMMON.FFIELD'
1319 include 'COMMON.SBRIDGE'
1320 include 'COMMON.HEADER'
1321 include 'COMMON.CONTROL'
1322 include 'COMMON.DBASE'
1323 include 'COMMON.THREAD'
1324 include 'COMMON.TIME1'
1325 include 'COMMON.SETUP'
1326 C Read bridging residues.
1327 read (inp,*) ns,(iss(i),i=1,ns)
1329 if(me.eq.king.or..not.out1file)
1330 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1331 C Check whether the specified bridging residues are cystines.
1333 if (itype(iss(i)).ne.1) then
1334 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1335 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1336 & ' can form a disulfide bridge?!!!'
1337 write (*,'(2a,i3,a)')
1338 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1339 & ' can form a disulfide bridge?!!!'
1341 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1346 C Read preformed bridges.
1348 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1349 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1352 C Check if the residues involved in bridges are in the specified list of
1353 C bridging residues.
1356 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1357 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1358 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1359 & ' contains residues present in other pairs.'
1360 write (*,'(a,i3,a)') 'Disulfide pair',i,
1361 & ' contains residues present in other pairs.'
1363 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1369 if (ihpb(i).eq.iss(j)) goto 10
1371 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1374 if (jhpb(i).eq.iss(j)) goto 20
1376 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1382 ihpb(i)=ihpb(i)+nres
1383 jhpb(i)=jhpb(i)+nres
1389 c----------------------------------------------------------------------------
1390 subroutine read_x(kanal,*)
1391 implicit real*8 (a-h,o-z)
1392 include 'DIMENSIONS'
1393 include 'COMMON.GEO'
1394 include 'COMMON.VAR'
1395 include 'COMMON.CHAIN'
1396 include 'COMMON.IOUNITS'
1397 include 'COMMON.CONTROL'
1398 include 'COMMON.LOCAL'
1399 include 'COMMON.INTERACT'
1400 c Read coordinates from input
1402 read(kanal,'(8f10.5)',end=10,err=10)
1403 & ((c(l,k),l=1,3),k=1,nres),
1404 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1407 c(j,2*nres)=c(j,nres)
1409 call int_from_cart1(.false.)
1412 dc(j,i)=c(j,i+1)-c(j,i)
1413 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1417 if (itype(i).ne.10) then
1419 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1420 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1428 c----------------------------------------------------------------------------
1429 subroutine read_threadbase
1430 implicit real*8 (a-h,o-z)
1431 include 'DIMENSIONS'
1432 include 'COMMON.IOUNITS'
1433 include 'COMMON.GEO'
1434 include 'COMMON.VAR'
1435 include 'COMMON.INTERACT'
1436 include 'COMMON.LOCAL'
1437 include 'COMMON.NAMES'
1438 include 'COMMON.CHAIN'
1439 include 'COMMON.FFIELD'
1440 include 'COMMON.SBRIDGE'
1441 include 'COMMON.HEADER'
1442 include 'COMMON.CONTROL'
1443 include 'COMMON.DBASE'
1444 include 'COMMON.THREAD'
1445 include 'COMMON.TIME1'
1446 C Read pattern database for threading.
1447 read (icbase,*) nseq
1449 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1450 & nres_base(2,i),nres_base(3,i)
1451 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1453 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1454 c & nres_base(2,i),nres_base(3,i)
1455 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1459 if (weidis.eq.0.0D0) weidis=0.1D0
1468 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1469 write (iout,'(a,i5)') 'nexcl: ',nexcl
1470 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1473 c------------------------------------------------------------------------------
1474 subroutine setup_var
1475 implicit real*8 (a-h,o-z)
1476 include 'DIMENSIONS'
1477 include 'COMMON.IOUNITS'
1478 include 'COMMON.GEO'
1479 include 'COMMON.VAR'
1480 include 'COMMON.INTERACT'
1481 include 'COMMON.LOCAL'
1482 include 'COMMON.NAMES'
1483 include 'COMMON.CHAIN'
1484 include 'COMMON.FFIELD'
1485 include 'COMMON.SBRIDGE'
1486 include 'COMMON.HEADER'
1487 include 'COMMON.CONTROL'
1488 include 'COMMON.DBASE'
1489 include 'COMMON.THREAD'
1490 include 'COMMON.TIME1'
1491 C Set up variable list.
1497 if (itype(i).ne.10) then
1499 ialph(i,1)=nvar+nside
1503 if (indphi.gt.0) then
1505 else if (indback.gt.0) then
1510 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1513 c----------------------------------------------------------------------------
1514 subroutine gen_dist_constr
1515 C Generate CA distance constraints.
1516 implicit real*8 (a-h,o-z)
1517 include 'DIMENSIONS'
1518 include 'COMMON.IOUNITS'
1519 include 'COMMON.GEO'
1520 include 'COMMON.VAR'
1521 include 'COMMON.INTERACT'
1522 include 'COMMON.LOCAL'
1523 include 'COMMON.NAMES'
1524 include 'COMMON.CHAIN'
1525 include 'COMMON.FFIELD'
1526 include 'COMMON.SBRIDGE'
1527 include 'COMMON.HEADER'
1528 include 'COMMON.CONTROL'
1529 include 'COMMON.DBASE'
1530 include 'COMMON.THREAD'
1531 include 'COMMON.TIME1'
1532 dimension itype_pdb(maxres)
1533 common /pizda/ itype_pdb
1535 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1536 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1537 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1539 do i=nstart_sup,nstart_sup+nsup-1
1540 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1541 cd & ' seq_pdb', restyp(itype_pdb(i))
1542 do j=i+2,nstart_sup+nsup-1
1544 ihpb(nhpb)=i+nstart_seq-nstart_sup
1545 jhpb(nhpb)=j+nstart_seq-nstart_sup
1547 dhpb(nhpb)=dist(i,j)
1550 cd write (iout,'(a)') 'Distance constraints:'
1555 cd if (ii.gt.nres) then
1560 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1561 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1562 cd & dhpb(i),forcon(i)
1566 c----------------------------------------------------------------------------
1568 implicit real*8 (a-h,o-z)
1569 include 'DIMENSIONS'
1570 include 'COMMON.MAP'
1571 include 'COMMON.IOUNITS'
1572 character*3 angid(4) /'THE','PHI','ALP','OME'/
1573 character*80 mapcard,ucase
1575 read (inp,'(a)') mapcard
1576 mapcard=ucase(mapcard)
1577 if (index(mapcard,'PHI').gt.0) then
1579 else if (index(mapcard,'THE').gt.0) then
1581 else if (index(mapcard,'ALP').gt.0) then
1583 else if (index(mapcard,'OME').gt.0) then
1586 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1587 stop 'Error - illegal variable spec in MAP card.'
1589 call readi (mapcard,'RES1',res1(imap),0)
1590 call readi (mapcard,'RES2',res2(imap),0)
1591 if (res1(imap).eq.0) then
1592 res1(imap)=res2(imap)
1593 else if (res2(imap).eq.0) then
1594 res2(imap)=res1(imap)
1596 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1598 & 'Error - illegal definition of variable group in MAP.'
1599 stop 'Error - illegal definition of variable group in MAP.'
1601 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1602 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1603 call readi(mapcard,'NSTEP',nstep(imap),0)
1604 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1606 & 'Illegal boundary and/or step size specification in MAP.'
1607 stop 'Illegal boundary and/or step size specification in MAP.'
1612 c----------------------------------------------------------------------------
1613 csa subroutine csaread
1614 csa implicit real*8 (a-h,o-z)
1615 csa include 'DIMENSIONS'
1616 csa include 'COMMON.IOUNITS'
1617 csa include 'COMMON.GEO'
1618 csa include 'COMMON.CSA'
1619 csa include 'COMMON.BANK'
1620 csa include 'COMMON.CONTROL'
1621 csa character*80 ucase
1622 csa character*620 mcmcard
1623 csa call card_concat(mcmcard)
1625 csa call readi(mcmcard,'NCONF',nconf,50)
1626 csa call readi(mcmcard,'NADD',nadd,0)
1627 csa call readi(mcmcard,'JSTART',jstart,1)
1628 csa call readi(mcmcard,'JEND',jend,1)
1629 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1630 csa call readi(mcmcard,'N0',n0,1)
1631 csa call readi(mcmcard,'N1',n1,6)
1632 csa call readi(mcmcard,'N2',n2,4)
1633 csa call readi(mcmcard,'N3',n3,0)
1634 csa call readi(mcmcard,'N4',n4,0)
1635 csa call readi(mcmcard,'N5',n5,0)
1636 csa call readi(mcmcard,'N6',n6,10)
1637 csa call readi(mcmcard,'N7',n7,0)
1638 csa call readi(mcmcard,'N8',n8,0)
1639 csa call readi(mcmcard,'N9',n9,0)
1640 csa call readi(mcmcard,'N14',n14,0)
1641 csa call readi(mcmcard,'N15',n15,0)
1642 csa call readi(mcmcard,'N16',n16,0)
1643 csa call readi(mcmcard,'N17',n17,0)
1644 csa call readi(mcmcard,'N18',n18,0)
1646 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1648 csa call readi(mcmcard,'NDIFF',ndiff,2)
1649 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1650 csa call readi(mcmcard,'IS1',is1,1)
1651 csa call readi(mcmcard,'IS2',is2,8)
1652 csa call readi(mcmcard,'NRAN0',nran0,4)
1653 csa call readi(mcmcard,'NRAN1',nran1,2)
1654 csa call readi(mcmcard,'IRR',irr,1)
1655 csa call readi(mcmcard,'NSEED',nseed,20)
1656 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1657 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1658 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1659 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1660 csa call readi(mcmcard,'ICMAX',icmax,3)
1661 csa call readi(mcmcard,'IRESTART',irestart,0)
1662 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1665 csa call reada(mcmcard,'DELE',dele,20.0d0)
1666 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1667 csa call readi(mcmcard,'IREF',iref,0)
1668 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1669 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1670 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1671 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1672 csa write (iout,*) "NCONF_IN",nconf_in
1675 c----------------------------------------------------------------------------
1676 cfmc subroutine mcmfread
1677 cfmc implicit real*8 (a-h,o-z)
1678 cfmc include 'DIMENSIONS'
1679 cfmc include 'COMMON.MCMF'
1680 cfmc include 'COMMON.IOUNITS'
1681 cfmc include 'COMMON.GEO'
1682 cfmc character*80 ucase
1683 cfmc character*620 mcmcard
1684 cfmc call card_concat(mcmcard)
1686 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1687 cfmc write(iout,*)'MAXRANT=',maxrant
1688 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1689 cfmc write(iout,*)'MAXFAM=',maxfam
1690 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1691 cfmc write(iout,*)'NNET1=',nnet1
1692 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1693 cfmc write(iout,*)'NNET2=',nnet2
1694 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1695 cfmc write(iout,*)'NNET3=',nnet3
1696 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1697 cfmc write(iout,*)'ILASTT=',ilastt
1698 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1699 cfmc write(iout,*)'MAXSTR=',maxstr
1700 cfmc maxstr_f=maxstr/maxfam
1701 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1702 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1703 cfmc write(iout,*)'NMCMF=',nmcmf
1704 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1705 cfmc write(iout,*)'IFOCUS=',ifocus
1706 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1707 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1708 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1709 cfmc write(iout,*)'INTPRT=',intprt
1710 cfmc call readi(mcmcard,'IPRT',iprt,100)
1711 cfmc write(iout,*)'IPRT=',iprt
1712 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1713 cfmc write(iout,*)'IMAXTR=',imaxtr
1714 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1715 cfmc write(iout,*)'MAXEVEN=',maxeven
1716 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1717 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1718 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1719 cfmc write(iout,*)'INIMIN=',inimin
1720 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1721 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1722 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1723 cfmc write(iout,*)'NTHREAD=',nthread
1724 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1725 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1726 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1727 cfmc write(iout,*)'MAXPERT=',maxpert
1728 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1729 cfmc write(iout,*)'IRMSD=',irmsd
1730 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1731 cfmc write(iout,*)'DENEMIN=',denemin
1732 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1733 cfmc write(iout,*)'RCUT1S=',rcut1s
1734 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1735 cfmc write(iout,*)'RCUT1E=',rcut1e
1736 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1737 cfmc write(iout,*)'RCUT2S=',rcut2s
1738 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1739 cfmc write(iout,*)'RCUT2E=',rcut2e
1740 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1741 cfmc write(iout,*)'DPERT1=',d_pert1
1742 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1743 cfmc write(iout,*)'DPERT1A=',d_pert1a
1744 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1745 cfmc write(iout,*)'DPERT2=',d_pert2
1746 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1747 cfmc write(iout,*)'DPERT2A=',d_pert2a
1748 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1749 cfmc write(iout,*)'DPERT2B=',d_pert2b
1750 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1751 cfmc write(iout,*)'DPERT2C=',d_pert2c
1752 cfmc d_pert1=deg2rad*d_pert1
1753 cfmc d_pert1a=deg2rad*d_pert1a
1754 cfmc d_pert2=deg2rad*d_pert2
1755 cfmc d_pert2a=deg2rad*d_pert2a
1756 cfmc d_pert2b=deg2rad*d_pert2b
1757 cfmc d_pert2c=deg2rad*d_pert2c
1758 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1759 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1760 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1761 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1762 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1763 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1764 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1765 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1766 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1767 cfmc write(iout,*)'RCUTINI=',rcutini
1768 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1769 cfmc write(iout,*)'GRAT=',grat
1770 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1771 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1775 c----------------------------------------------------------------------------
1777 implicit real*8 (a-h,o-z)
1778 include 'DIMENSIONS'
1779 include 'COMMON.MCM'
1780 include 'COMMON.MCE'
1781 include 'COMMON.IOUNITS'
1783 character*320 mcmcard
1784 call card_concat(mcmcard)
1785 call readi(mcmcard,'MAXACC',maxacc,100)
1786 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1787 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1788 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1789 call readi(mcmcard,'MAXREPM',maxrepm,200)
1790 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1791 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1792 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1793 call reada(mcmcard,'E_UP',e_up,5.0D0)
1794 call reada(mcmcard,'DELTE',delte,0.1D0)
1795 call readi(mcmcard,'NSWEEP',nsweep,5)
1796 call readi(mcmcard,'NSTEPH',nsteph,0)
1797 call readi(mcmcard,'NSTEPC',nstepc,0)
1798 call reada(mcmcard,'TMIN',tmin,298.0D0)
1799 call reada(mcmcard,'TMAX',tmax,298.0D0)
1800 call readi(mcmcard,'NWINDOW',nwindow,0)
1801 call readi(mcmcard,'PRINT_MC',print_mc,0)
1802 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1803 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1804 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1805 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1806 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1807 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1808 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1809 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1810 if (nwindow.gt.0) then
1811 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1813 winlen(i)=winend(i)-winstart(i)+1
1816 if (tmax.lt.tmin) tmax=tmin
1817 if (tmax.eq.tmin) then
1821 if (nstepc.gt.0 .and. nsteph.gt.0) then
1822 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1823 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1825 C Probabilities of different move types
1826 sumpro_type(0)=0.0D0
1827 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1828 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1829 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1830 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1831 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1832 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1833 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1835 print *,'i',i,' sumprotype',sumpro_type(i)
1836 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1837 print *,'i',i,' sumprotype',sumpro_type(i)
1841 c----------------------------------------------------------------------------
1842 subroutine read_minim
1843 implicit real*8 (a-h,o-z)
1844 include 'DIMENSIONS'
1845 include 'COMMON.MINIM'
1846 include 'COMMON.IOUNITS'
1848 character*320 minimcard
1849 call card_concat(minimcard)
1850 call readi(minimcard,'MAXMIN',maxmin,2000)
1851 call readi(minimcard,'MAXFUN',maxfun,5000)
1852 call readi(minimcard,'MINMIN',minmin,maxmin)
1853 call readi(minimcard,'MINFUN',minfun,maxmin)
1854 call reada(minimcard,'TOLF',tolf,1.0D-2)
1855 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1856 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1857 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1858 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1859 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1860 & 'Options in energy minimization:'
1861 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1862 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1863 & 'MinMin:',MinMin,' MinFun:',MinFun,
1864 & ' TolF:',TolF,' RTolF:',RTolF
1867 c----------------------------------------------------------------------------
1868 subroutine read_angles(kanal,*)
1869 implicit real*8 (a-h,o-z)
1870 include 'DIMENSIONS'
1871 include 'COMMON.GEO'
1872 include 'COMMON.VAR'
1873 include 'COMMON.CHAIN'
1874 include 'COMMON.IOUNITS'
1875 include 'COMMON.CONTROL'
1876 c Read angles from input
1878 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1879 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1880 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1881 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1884 c 9/7/01 avoid 180 deg valence angle
1885 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1887 theta(i)=deg2rad*theta(i)
1888 phi(i)=deg2rad*phi(i)
1889 alph(i)=deg2rad*alph(i)
1890 omeg(i)=deg2rad*omeg(i)
1895 c----------------------------------------------------------------------------
1896 subroutine reada(rekord,lancuch,wartosc,default)
1898 character*(*) rekord,lancuch
1899 double precision wartosc,default
1902 iread=index(rekord,lancuch)
1903 if (iread.eq.0) then
1907 iread=iread+ilen(lancuch)+1
1908 read (rekord(iread:),*,err=10,end=10) wartosc
1913 c----------------------------------------------------------------------------
1914 subroutine readi(rekord,lancuch,wartosc,default)
1916 character*(*) rekord,lancuch
1917 integer wartosc,default
1920 iread=index(rekord,lancuch)
1921 if (iread.eq.0) then
1925 iread=iread+ilen(lancuch)+1
1926 read (rekord(iread:),*,err=10,end=10) wartosc
1931 c----------------------------------------------------------------------------
1932 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1935 integer tablica(dim),default
1936 character*(*) rekord,lancuch
1943 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1944 if (iread.eq.0) return
1945 iread=iread+ilen(lancuch)+1
1946 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1949 c----------------------------------------------------------------------------
1950 subroutine multreada(rekord,lancuch,tablica,dim,default)
1953 double precision tablica(dim),default
1954 character*(*) rekord,lancuch
1961 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1962 if (iread.eq.0) return
1963 iread=iread+ilen(lancuch)+1
1964 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1967 c----------------------------------------------------------------------------
1968 subroutine openunits
1969 implicit real*8 (a-h,o-z)
1970 include 'DIMENSIONS'
1973 character*16 form,nodename
1976 include 'COMMON.SETUP'
1977 include 'COMMON.IOUNITS'
1979 include 'COMMON.CONTROL'
1980 integer lenpre,lenpot,ilen,lentmp
1982 character*3 out1file_text,ucase
1985 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1986 call getenv_loc("PREFIX",prefix)
1988 call getenv_loc("POT",pot)
1989 call getenv_loc("DIRTMP",tmpdir)
1990 call getenv_loc("CURDIR",curdir)
1991 call getenv_loc("OUT1FILE",out1file_text)
1992 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1993 out1file_text=ucase(out1file_text)
1994 if (out1file_text(1:1).eq."Y") then
1997 out1file=fg_rank.gt.0
2002 if (lentmp.gt.0) then
2003 write (*,'(80(1h!))')
2004 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2005 write (*,'(80(1h!))')
2006 write (*,*)"All output files will be on node /tmp directory."
2008 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2009 if (me.eq.king) then
2010 write (*,*) "The master node is ",nodename
2011 else if (fg_rank.eq.0) then
2012 write (*,*) "I am the CG slave node ",nodename
2014 write (*,*) "I am the FG slave node ",nodename
2017 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2018 lenpre = lentmp+lenpre+1
2020 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2021 C Get the names and open the input files
2022 #if defined(WINIFL) || defined(WINPGI)
2023 open(1,file=pref_orig(:ilen(pref_orig))//
2024 & '.inp',status='old',readonly,shared)
2025 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2026 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2027 C Get parameter filenames and open the parameter files.
2028 call getenv_loc('BONDPAR',bondname)
2029 open (ibond,file=bondname,status='old',readonly,shared)
2030 call getenv_loc('THETPAR',thetname)
2031 open (ithep,file=thetname,status='old',readonly,shared)
2033 call getenv_loc('THETPARPDB',thetname_pdb)
2034 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2036 call getenv_loc('ROTPAR',rotname)
2037 open (irotam,file=rotname,status='old',readonly,shared)
2039 call getenv_loc('ROTPARPDB',rotname_pdb)
2040 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2042 call getenv_loc('TORPAR',torname)
2043 open (itorp,file=torname,status='old',readonly,shared)
2044 call getenv_loc('TORDPAR',tordname)
2045 open (itordp,file=tordname,status='old',readonly,shared)
2046 call getenv_loc('FOURIER',fouriername)
2047 open (ifourier,file=fouriername,status='old',readonly,shared)
2048 call getenv_loc('ELEPAR',elename)
2049 open (ielep,file=elename,status='old',readonly,shared)
2050 call getenv_loc('SIDEPAR',sidename)
2051 open (isidep,file=sidename,status='old',readonly,shared)
2052 #elif (defined CRAY) || (defined AIX)
2053 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2055 c print *,"Processor",myrank," opened file 1"
2056 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2057 c print *,"Processor",myrank," opened file 9"
2058 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2059 C Get parameter filenames and open the parameter files.
2060 call getenv_loc('BONDPAR',bondname)
2061 open (ibond,file=bondname,status='old',action='read')
2062 c print *,"Processor",myrank," opened file IBOND"
2063 call getenv_loc('THETPAR',thetname)
2064 open (ithep,file=thetname,status='old',action='read')
2065 c print *,"Processor",myrank," opened file ITHEP"
2067 call getenv_loc('THETPARPDB',thetname_pdb)
2068 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2070 call getenv_loc('ROTPAR',rotname)
2071 open (irotam,file=rotname,status='old',action='read')
2072 c print *,"Processor",myrank," opened file IROTAM"
2074 call getenv_loc('ROTPARPDB',rotname_pdb)
2075 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2077 call getenv_loc('TORPAR',torname)
2078 open (itorp,file=torname,status='old',action='read')
2079 c print *,"Processor",myrank," opened file ITORP"
2080 call getenv_loc('TORDPAR',tordname)
2081 open (itordp,file=tordname,status='old',action='read')
2082 c print *,"Processor",myrank," opened file ITORDP"
2083 call getenv_loc('SCCORPAR',sccorname)
2084 open (isccor,file=sccorname,status='old',action='read')
2085 c print *,"Processor",myrank," opened file ISCCOR"
2086 call getenv_loc('FOURIER',fouriername)
2087 open (ifourier,file=fouriername,status='old',action='read')
2088 c print *,"Processor",myrank," opened file IFOURIER"
2089 call getenv_loc('ELEPAR',elename)
2090 open (ielep,file=elename,status='old',action='read')
2091 c print *,"Processor",myrank," opened file IELEP"
2092 call getenv_loc('SIDEPAR',sidename)
2093 open (isidep,file=sidename,status='old',action='read')
2094 c print *,"Processor",myrank," opened file ISIDEP"
2095 c print *,"Processor",myrank," opened parameter files"
2097 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2098 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2099 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2100 C Get parameter filenames and open the parameter files.
2101 call getenv_loc('BONDPAR',bondname)
2102 open (ibond,file=bondname,status='old')
2103 call getenv_loc('THETPAR',thetname)
2104 open (ithep,file=thetname,status='old')
2106 call getenv_loc('THETPARPDB',thetname_pdb)
2107 open (ithep_pdb,file=thetname_pdb,status='old')
2109 call getenv_loc('ROTPAR',rotname)
2110 open (irotam,file=rotname,status='old')
2112 call getenv_loc('ROTPARPDB',rotname_pdb)
2113 open (irotam_pdb,file=rotname_pdb,status='old')
2115 call getenv_loc('TORPAR',torname)
2116 open (itorp,file=torname,status='old')
2117 call getenv_loc('TORDPAR',tordname)
2118 open (itordp,file=tordname,status='old')
2119 call getenv_loc('SCCORPAR',sccorname)
2120 open (isccor,file=sccorname,status='old')
2121 call getenv_loc('FOURIER',fouriername)
2122 open (ifourier,file=fouriername,status='old')
2123 call getenv_loc('ELEPAR',elename)
2124 open (ielep,file=elename,status='old')
2125 call getenv_loc('SIDEPAR',sidename)
2126 open (isidep,file=sidename,status='old')
2128 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2130 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2131 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2132 C Get parameter filenames and open the parameter files.
2133 call getenv_loc('BONDPAR',bondname)
2134 open (ibond,file=bondname,status='old',action='read')
2135 call getenv_loc('THETPAR',thetname)
2136 open (ithep,file=thetname,status='old',action='read')
2138 call getenv_loc('THETPARPDB',thetname_pdb)
2139 print *,"thetname_pdb ",thetname_pdb
2140 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2141 print *,ithep_pdb," opened"
2143 call getenv_loc('ROTPAR',rotname)
2144 open (irotam,file=rotname,status='old',action='read')
2146 call getenv_loc('ROTPARPDB',rotname_pdb)
2147 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2149 call getenv_loc('TORPAR',torname)
2150 open (itorp,file=torname,status='old',action='read')
2151 call getenv_loc('TORDPAR',tordname)
2152 open (itordp,file=tordname,status='old',action='read')
2153 call getenv_loc('SCCORPAR',sccorname)
2154 open (isccor,file=sccorname,status='old',action='read')
2155 call getenv_loc('FOURIER',fouriername)
2156 open (ifourier,file=fouriername,status='old',action='read')
2157 call getenv_loc('ELEPAR',elename)
2158 open (ielep,file=elename,status='old',action='read')
2159 call getenv_loc('SIDEPAR',sidename)
2160 open (isidep,file=sidename,status='old',action='read')
2164 C 8/9/01 In the newest version SCp interaction constants are read from a file
2165 C Use -DOLDSCP to use hard-coded constants instead.
2167 call getenv_loc('SCPPAR',scpname)
2168 #if defined(WINIFL) || defined(WINPGI)
2169 open (iscpp,file=scpname,status='old',readonly,shared)
2170 #elif (defined CRAY) || (defined AIX)
2171 open (iscpp,file=scpname,status='old',action='read')
2173 open (iscpp,file=scpname,status='old')
2175 open (iscpp,file=scpname,status='old',action='read')
2178 call getenv_loc('PATTERN',patname)
2179 #if defined(WINIFL) || defined(WINPGI)
2180 open (icbase,file=patname,status='old',readonly,shared)
2181 #elif (defined CRAY) || (defined AIX)
2182 open (icbase,file=patname,status='old',action='read')
2184 open (icbase,file=patname,status='old')
2186 open (icbase,file=patname,status='old',action='read')
2189 C Open output file only for CG processes
2190 c print *,"Processor",myrank," fg_rank",fg_rank
2191 if (fg_rank.eq.0) then
2193 if (nodes.eq.1) then
2196 npos = dlog10(dfloat(nodes-1))+1
2198 if (npos.lt.3) npos=3
2199 write (liczba,'(i1)') npos
2200 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2202 write (liczba,form) me
2203 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2204 & liczba(:ilen(liczba))
2205 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2207 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2209 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2210 & liczba(:ilen(liczba))//'.mol2'
2211 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2212 & liczba(:ilen(liczba))//'.stat'
2214 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2215 & //liczba(:ilen(liczba))//'.stat')
2216 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2219 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2220 & liczba(:ilen(liczba))//'.const'
2225 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2226 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2227 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2228 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2229 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2231 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2233 rest2name=prefix(:ilen(prefix))//'.rst'
2235 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2238 #if defined(AIX) || defined(PGI)
2239 if (me.eq.king .or. .not. out1file)
2240 & open(iout,file=outname,status='unknown')
2243 if (fg_rank.gt.0) then
2244 write (liczba,'(i3.3)') myrank/nfgtasks
2245 write (ll,'(bz,i3.3)') fg_rank
2246 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2252 open(igeom,file=intname,status='unknown',position='append')
2253 open(ipdb,file=pdbname,status='unknown')
2254 open(imol2,file=mol2name,status='unknown')
2255 open(istat,file=statname,status='unknown',position='append')
2257 c1out open(iout,file=outname,status='unknown')
2260 if (me.eq.king .or. .not.out1file)
2261 & open(iout,file=outname,status='unknown')
2264 if (fg_rank.gt.0) then
2265 print "Processor",fg_rank," opening output file"
2266 write (liczba,'(i3.3)') myrank/nfgtasks
2267 write (ll,'(bz,i3.3)') fg_rank
2268 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2274 open(igeom,file=intname,status='unknown',access='append')
2275 open(ipdb,file=pdbname,status='unknown')
2276 open(imol2,file=mol2name,status='unknown')
2277 open(istat,file=statname,status='unknown',access='append')
2279 c1out open(iout,file=outname,status='unknown')
2282 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2283 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2284 csa csa_history=prefix(:lenpre)//'.CSA.history'
2285 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2286 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2287 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2288 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2289 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2290 csa csa_int=prefix(:lenpre)//'.int'
2291 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2292 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2293 csa csa_in=prefix(:lenpre)//'.CSA.in'
2294 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2297 write (iout,'(80(1h-))')
2298 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2299 write (iout,'(80(1h-))')
2300 write (iout,*) "Input file : ",
2301 & pref_orig(:ilen(pref_orig))//'.inp'
2302 write (iout,*) "Output file : ",
2303 & outname(:ilen(outname))
2305 write (iout,*) "Sidechain potential file : ",
2306 & sidename(:ilen(sidename))
2308 write (iout,*) "SCp potential file : ",
2309 & scpname(:ilen(scpname))
2311 write (iout,*) "Electrostatic potential file : ",
2312 & elename(:ilen(elename))
2313 write (iout,*) "Cumulant coefficient file : ",
2314 & fouriername(:ilen(fouriername))
2315 write (iout,*) "Torsional parameter file : ",
2316 & torname(:ilen(torname))
2317 write (iout,*) "Double torsional parameter file : ",
2318 & tordname(:ilen(tordname))
2319 write (iout,*) "SCCOR parameter file : ",
2320 & sccorname(:ilen(sccorname))
2321 write (iout,*) "Bond & inertia constant file : ",
2322 & bondname(:ilen(bondname))
2323 write (iout,*) "Bending parameter file : ",
2324 & thetname(:ilen(thetname))
2325 write (iout,*) "Rotamer parameter file : ",
2326 & rotname(:ilen(rotname))
2327 write (iout,*) "Threading database : ",
2328 & patname(:ilen(patname))
2330 &write (iout,*)" DIRTMP : ",
2332 write (iout,'(80(1h-))')
2336 c----------------------------------------------------------------------------
2337 subroutine card_concat(card)
2338 implicit real*8 (a-h,o-z)
2339 include 'DIMENSIONS'
2340 include 'COMMON.IOUNITS'
2342 character*80 karta,ucase
2344 read (inp,'(a)') karta
2347 do while (karta(80:80).eq.'&')
2348 card=card(:ilen(card)+1)//karta(:79)
2349 read (inp,'(a)') karta
2352 card=card(:ilen(card)+1)//karta
2355 c----------------------------------------------------------------------------------
2357 implicit real*8 (a-h,o-z)
2358 include 'DIMENSIONS'
2359 include 'COMMON.CHAIN'
2360 include 'COMMON.IOUNITS'
2362 open(irest2,file=rest2name,status='unknown')
2363 read(irest2,*) totT,EK,potE,totE,t_bath
2365 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2368 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2371 read (irest2,*) iset
2376 c---------------------------------------------------------------------------------
2377 subroutine read_fragments
2378 implicit real*8 (a-h,o-z)
2379 include 'DIMENSIONS'
2383 include 'COMMON.SETUP'
2384 include 'COMMON.CHAIN'
2385 include 'COMMON.IOUNITS'
2387 include 'COMMON.CONTROL'
2388 read(inp,*) nset,nfrag,npair,nfrag_back
2389 if(me.eq.king.or..not.out1file)
2390 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2391 & " nfrag_back",nfrag_back
2393 read(inp,*) mset(iset)
2395 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2397 if(me.eq.king.or..not.out1file)
2398 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2399 & ifrag(2,i,iset), qinfrag(i,iset)
2402 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2404 if(me.eq.king.or..not.out1file)
2405 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2406 & ipair(2,i,iset), qinpair(i,iset)
2409 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2410 & wfrag_back(3,i,iset),
2411 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2412 if(me.eq.king.or..not.out1file)
2413 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2414 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2419 c-------------------------------------------------------------------------------
2420 subroutine read_dist_constr
2421 implicit real*8 (a-h,o-z)
2422 include 'DIMENSIONS'
2426 include 'COMMON.SETUP'
2427 include 'COMMON.CONTROL'
2428 include 'COMMON.CHAIN'
2429 include 'COMMON.IOUNITS'
2430 include 'COMMON.SBRIDGE'
2431 integer ifrag_(2,100),ipair_(2,100)
2432 double precision wfrag_(100),wpair_(100)
2433 character*500 controlcard
2434 c write (iout,*) "Calling read_dist_constr"
2435 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2437 call card_concat(controlcard)
2438 call readi(controlcard,"NFRAG",nfrag_,0)
2439 call readi(controlcard,"NPAIR",npair_,0)
2440 call readi(controlcard,"NDIST",ndist_,0)
2441 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2442 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2443 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2444 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2445 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2446 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2447 c write (iout,*) "IFRAG"
2449 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2451 c write (iout,*) "IPAIR"
2453 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2455 if (.not.refstr .and. nfrag.gt.0) then
2457 & "ERROR: no reference structure to compute distance restraints"
2459 & "Restraints must be specified explicitly (NDIST=number)"
2462 if (nfrag.lt.2 .and. npair.gt.0) then
2463 write (iout,*) "ERROR: Less than 2 fragments specified",
2464 & " but distance restraints between pairs requested"
2469 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2470 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2471 & ifrag_(2,i)=nstart_sup+nsup-1
2472 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2474 if (wfrag_(i).gt.0.0d0) then
2475 do j=ifrag_(1,i),ifrag_(2,i)-1
2476 do k=j+1,ifrag_(2,i)
2477 write (iout,*) "j",j," k",k
2479 if (constr_dist.eq.1) then
2484 forcon(nhpb)=wfrag_(i)
2485 else if (constr_dist.eq.2) then
2486 if (ddjk.le.dist_cut) then
2491 forcon(nhpb)=wfrag_(i)
2498 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2501 if (.not.out1file .or. me.eq.king)
2502 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2503 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2505 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2506 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2513 if (wpair_(i).gt.0.0d0) then
2521 do j=ifrag_(1,ii),ifrag_(2,ii)
2522 do k=ifrag_(1,jj),ifrag_(2,jj)
2526 forcon(nhpb)=wpair_(i)
2527 dhpb(nhpb)=dist(j,k)
2529 if (.not.out1file .or. me.eq.king)
2530 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2531 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2533 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2534 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2541 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2542 & ibecarb(i),forcon(nhpb+1)
2543 if (forcon(nhpb+1).gt.0.0d0) then
2545 if (ibecarb(i).gt.0) then
2546 ihpb(i)=ihpb(i)+nres
2547 jhpb(i)=jhpb(i)+nres
2549 if (dhpb(nhpb).eq.0.0d0)
2550 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2554 if (.not.out1file .or. me.eq.king) then
2557 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2558 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2566 c-------------------------------------------------------------------------------
2568 subroutine flush(iu)
2573 subroutine flush(iu)
2578 c------------------------------------------------------------------------------
2579 subroutine copy_to_tmp(source)
2580 include "DIMENSIONS"
2581 include "COMMON.IOUNITS"
2582 character*(*) source
2583 character* 256 tmpfile
2587 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2588 inquire(file=tmpfile,exist=ex)
2590 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2591 & " to temporary directory..."
2592 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2593 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2597 c------------------------------------------------------------------------------
2598 subroutine move_from_tmp(source)
2599 include "DIMENSIONS"
2600 include "COMMON.IOUNITS"
2601 character*(*) source
2604 write (*,*) "Moving ",source(:ilen(source)),
2605 & " from temporary directory to working directory"
2606 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2607 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2610 c------------------------------------------------------------------------------
2611 subroutine random_init(seed)
2613 C Initialize random number generator
2615 implicit real*8 (a-h,o-z)
2616 include 'DIMENSIONS'
2622 logical OKRandom, prng_restart
2624 integer iseed_array(4)
2626 include 'COMMON.IOUNITS'
2627 include 'COMMON.TIME1'
2628 include 'COMMON.THREAD'
2629 include 'COMMON.SBRIDGE'
2630 include 'COMMON.CONTROL'
2631 include 'COMMON.MCM'
2632 include 'COMMON.MAP'
2633 include 'COMMON.HEADER'
2634 csa include 'COMMON.CSA'
2635 include 'COMMON.CHAIN'
2636 include 'COMMON.MUCA'
2638 include 'COMMON.FFIELD'
2639 include 'COMMON.SETUP'
2640 iseed=-dint(dabs(seed))
2641 if (iseed.eq.0) then
2642 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2643 & 'Random seed undefined. The program will stop.'
2644 write (*,'(/80(1h*)/20x,a/80(1h*))')
2645 & 'Random seed undefined. The program will stop.'
2647 call mpi_finalize(mpi_comm_world,ierr)
2649 stop 'Bad random seed.'
2652 if (fg_rank.eq.0) then
2656 if(me.eq.king .or. .not. out1file)
2657 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2658 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2659 OKRandom = prng_restart(me,iseedi8)
2662 tmp=65536.0d0**(4-i)
2663 iseed_array(i) = dint(seed/tmp)
2664 seed=seed-iseed_array(i)*tmp
2666 if(me.eq.king .or. .not. out1file)
2667 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2668 & (iseed_array(i),i=1,4)
2669 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2670 & (iseed_array(i),i=1,4)
2671 OKRandom = prng_restart(me,iseed_array)
2674 r1=ran_number(0.0D0,1.0D0)
2675 if(me.eq.king .or. .not. out1file)
2676 & write (iout,*) 'ran_num',r1
2677 if (r1.lt.0.0d0) OKRandom=.false.
2679 if (.not.OKRandom) then
2680 write (iout,*) 'PRNG IS NOT WORKING!!!'
2681 print *,'PRNG IS NOT WORKING!!!'
2684 call mpi_abort(mpi_comm_world,error_msg,ierr)
2687 write (iout,*) 'too many processors for parallel prng'
2688 write (*,*) 'too many processors for parallel prng'
2696 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)