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 C Following 2 lines for diagnostics; comment out if not needed
921 write (iout,*) "Before sideadd"
923 if(me.eq.king.or..not.out1file)
924 & write(iout,*)'Adding sidechains'
931 do while (fail.and.nsi.le.maxsi)
932 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
935 if(fail) write(iout,*)'Adding sidechain failed for res ',
936 & i,' after ',nsi,' trials'
939 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
942 C Following 2 lines for diagnostics; comment out if not needed
943 write (iout,*) "After sideadd"
946 if (indpdb.eq.0) then
947 C Read sequence if not taken from the pdb file.
949 c print *,'nres=',nres
950 if (iscode.gt.0) then
951 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
953 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
955 C Convert sequence to numeric code
957 itype(i)=rescode(i,sequence(i),iscode)
959 C Assign initial virtual bond lengths
965 vbld(i+nres)=dsc(itype(i))
966 vbld_inv(i+nres)=dsc_inv(itype(i))
967 c write (iout,*) "i",i," itype",itype(i),
968 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
972 c print '(20i4)',(itype(i),i=1,nres)
975 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
977 if (itype(i).eq.21) then
981 else if (itype(i+1).ne.20) then
983 else if (itype(i).ne.20) then
990 if(me.eq.king.or..not.out1file)then
991 write (iout,*) "ITEL"
993 write (iout,*) i,itype(i),itel(i)
995 print *,'Call Read_Bridge.'
998 C 8/13/98 Set limits to generating the dihedral angles
1003 read (inp,*) ndih_constr
1004 if (ndih_constr.gt.0) then
1006 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1007 if(me.eq.king.or..not.out1file)then
1009 & 'There are',ndih_constr,' constraints on phi angles.'
1011 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1015 phi0(i)=deg2rad*phi0(i)
1016 drange(i)=deg2rad*drange(i)
1018 if(me.eq.king.or..not.out1file)
1019 & write (iout,*) 'FTORS',ftors
1022 phibound(1,ii) = phi0(i)-drange(i)
1023 phibound(2,ii) = phi0(i)+drange(i)
1028 if (me.eq.king) then
1030 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1032 write (iout,'(a3,i5,2f10.1)')
1033 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1039 cd print *,'NNT=',NNT,' NCT=',NCT
1040 if (itype(1).eq.21) nnt=2
1041 if (itype(nres).eq.21) nct=nct-1
1043 if(me.eq.king.or..not.out1file)
1044 & write (iout,'(a,i3)') 'nsup=',nsup
1046 if (nsup.le.(nct-nnt+1)) then
1047 do i=0,nct-nnt+1-nsup
1048 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1054 & 'Error - sequences to be superposed do not match.'
1057 do i=0,nsup-(nct-nnt+1)
1058 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1060 nstart_sup=nstart_sup+i
1066 & 'Error - sequences to be superposed do not match.'
1069 if (nsup.eq.0) nsup=nct-nnt
1070 if (nstart_sup.eq.0) nstart_sup=nnt
1071 if (nstart_seq.eq.0) nstart_seq=nnt
1072 if(me.eq.king.or..not.out1file)
1073 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1074 & ' nstart_seq=',nstart_seq
1076 c--- Zscore rms -------
1077 if (nz_start.eq.0) nz_start=nnt
1078 if (nz_end.eq.0 .and. nsup.gt.0) then
1080 else if (nz_end.eq.0) then
1083 if(me.eq.king.or..not.out1file)then
1084 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1085 write (iout,*) 'IZ_SC=',iz_sc
1087 c----------------------
1090 if (.not.pdbref) then
1091 call read_angles(inp,*38)
1093 38 write (iout,'(a)') 'Error reading reference structure.'
1095 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1096 stop 'Error reading reference structure'
1100 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1109 call contact(.true.,ncont_ref,icont_ref,co)
1111 if(me.eq.king.or..not.out1file)
1112 & write (iout,*) 'Contact order:',co
1114 if(me.eq.king.or..not.out1file)
1115 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1118 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1120 if(me.eq.king.or..not.out1file)
1121 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1122 & icont_ref(1,i),' ',
1123 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1127 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1128 if (constr_dist.gt.0) then
1129 call read_dist_constr
1131 if (nhpb.gt.0) call hpb_partition
1132 c write (iout,*) "After read_dist_constr nhpb",nhpb
1134 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1135 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1136 & modecalc.ne.10) then
1137 C If input structure hasn't been supplied from the PDB file read or generate
1139 if (iranconf.eq.0 .and. .not. extconf) then
1140 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1141 & write (iout,'(a)') 'Initial geometry will be read in.'
1143 read(inp,'(8f10.5)',end=36,err=36)
1144 & ((c(l,k),l=1,3),k=1,nres),
1145 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1146 call int_from_cart1(.false.)
1149 dc(j,i)=c(j,i+1)-c(j,i)
1150 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1154 if (itype(i).ne.10) then
1156 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1157 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1163 call read_angles(inp,*36)
1166 36 write (iout,'(a)') 'Error reading angle file.'
1168 call mpi_finalize( MPI_COMM_WORLD,IERR )
1170 stop 'Error reading angle file.'
1172 else if (extconf) then
1173 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1174 & write (iout,'(a)') 'Extended chain initial geometry.'
1176 theta(i)=90d0*deg2rad
1179 phi(i)=180d0*deg2rad
1182 alph(i)=110d0*deg2rad
1185 omeg(i)=-120d0*deg2rad
1188 if(me.eq.king.or..not.out1file)
1189 & write (iout,'(a)') 'Random-generated initial geometry.'
1193 if (me.eq.king .or. fg_rank.eq.0 .and. (
1194 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1198 call gen_rand_conf(itmp,*30)
1200 30 write (iout,*) 'Failed to generate random conformation',
1201 & ', itrial=',itrial
1202 write (*,*) 'Processor:',me,
1203 & ' Failed to generate random conformation',
1213 write (iout,'(a,i3,a)') 'Processor:',me,
1214 & ' error in generating random conformation.'
1215 write (*,'(a,i3,a)') 'Processor:',me,
1216 & ' error in generating random conformation.'
1219 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1226 elseif (modecalc.eq.4) then
1227 read (inp,'(a)') intinname
1228 open (intin,file=intinname,status='old',err=333)
1229 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1230 & write (iout,'(a)') 'intinname',intinname
1231 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1233 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1235 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1237 stop 'Error opening angle file.'
1241 C Generate distance constraints, if the PDB structure is to be regularized.
1242 if (nthread.gt.0) then
1243 call read_threadbase
1246 if (me.eq.king .or. .not. out1file)
1248 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1249 write (iout,'(/a,i3,a)')
1250 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1251 write (iout,'(20i4)') (iss(i),i=1,ns)
1253 write(iout,*)"Running with dynamic disulfide-bond formation"
1257 forcon(i-nss)=forcon(i)
1264 dyn_ss_mask(iss(i))=.true.
1267 write (iout,'(/a/)') 'Pre-formed links are:'
1273 if (me.eq.king.or..not.out1file)
1274 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1275 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1281 if (i2ndstr.gt.0) call secstrp2dihc
1282 c call geom_to_var(nvar,x)
1283 c call etotal(energia(0))
1284 c call enerprint(energia(0))
1285 c call briefout(0,etot)
1287 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1288 cd write (iout,'(a)') 'Variable list:'
1289 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1291 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1292 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1293 & 'Processor',myrank,': end reading molecular data.'
1297 c--------------------------------------------------------------------------
1298 logical function seq_comp(itypea,itypeb,length)
1300 integer length,itypea(length),itypeb(length)
1303 if (itypea(i).ne.itypeb(i)) then
1311 c-----------------------------------------------------------------------------
1312 subroutine read_bridge
1313 C Read information about disulfide bridges.
1314 implicit real*8 (a-h,o-z)
1315 include 'DIMENSIONS'
1319 include 'COMMON.IOUNITS'
1320 include 'COMMON.GEO'
1321 include 'COMMON.VAR'
1322 include 'COMMON.INTERACT'
1323 include 'COMMON.LOCAL'
1324 include 'COMMON.NAMES'
1325 include 'COMMON.CHAIN'
1326 include 'COMMON.FFIELD'
1327 include 'COMMON.SBRIDGE'
1328 include 'COMMON.HEADER'
1329 include 'COMMON.CONTROL'
1330 include 'COMMON.DBASE'
1331 include 'COMMON.THREAD'
1332 include 'COMMON.TIME1'
1333 include 'COMMON.SETUP'
1334 C Read bridging residues.
1335 read (inp,*) ns,(iss(i),i=1,ns)
1337 if(me.eq.king.or..not.out1file)
1338 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1339 C Check whether the specified bridging residues are cystines.
1341 if (itype(iss(i)).ne.1) then
1342 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1343 & 'Do you REALLY think that the residue ',
1344 & restyp(itype(iss(i))),i,
1345 & ' can form a disulfide bridge?!!!'
1346 write (*,'(2a,i3,a)')
1347 & 'Do you REALLY think that the residue ',
1348 & restyp(itype(iss(i))),i,
1349 & ' can form a disulfide bridge?!!!'
1351 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1356 C Read preformed bridges.
1358 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1359 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1362 C Check if the residues involved in bridges are in the specified list of
1363 C bridging residues.
1366 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1367 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1368 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1369 & ' contains residues present in other pairs.'
1370 write (*,'(a,i3,a)') 'Disulfide pair',i,
1371 & ' contains residues present in other pairs.'
1373 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1379 if (ihpb(i).eq.iss(j)) goto 10
1381 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1384 if (jhpb(i).eq.iss(j)) goto 20
1386 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1392 ihpb(i)=ihpb(i)+nres
1393 jhpb(i)=jhpb(i)+nres
1399 c----------------------------------------------------------------------------
1400 subroutine read_x(kanal,*)
1401 implicit real*8 (a-h,o-z)
1402 include 'DIMENSIONS'
1403 include 'COMMON.GEO'
1404 include 'COMMON.VAR'
1405 include 'COMMON.CHAIN'
1406 include 'COMMON.IOUNITS'
1407 include 'COMMON.CONTROL'
1408 include 'COMMON.LOCAL'
1409 include 'COMMON.INTERACT'
1410 c Read coordinates from input
1412 read(kanal,'(8f10.5)',end=10,err=10)
1413 & ((c(l,k),l=1,3),k=1,nres),
1414 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1417 c(j,2*nres)=c(j,nres)
1419 call int_from_cart1(.false.)
1422 dc(j,i)=c(j,i+1)-c(j,i)
1423 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1427 if (itype(i).ne.10) then
1429 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1430 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1438 c----------------------------------------------------------------------------
1439 subroutine read_threadbase
1440 implicit real*8 (a-h,o-z)
1441 include 'DIMENSIONS'
1442 include 'COMMON.IOUNITS'
1443 include 'COMMON.GEO'
1444 include 'COMMON.VAR'
1445 include 'COMMON.INTERACT'
1446 include 'COMMON.LOCAL'
1447 include 'COMMON.NAMES'
1448 include 'COMMON.CHAIN'
1449 include 'COMMON.FFIELD'
1450 include 'COMMON.SBRIDGE'
1451 include 'COMMON.HEADER'
1452 include 'COMMON.CONTROL'
1453 include 'COMMON.DBASE'
1454 include 'COMMON.THREAD'
1455 include 'COMMON.TIME1'
1456 C Read pattern database for threading.
1457 read (icbase,*) nseq
1459 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1460 & nres_base(2,i),nres_base(3,i)
1461 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1463 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1464 c & nres_base(2,i),nres_base(3,i)
1465 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1469 if (weidis.eq.0.0D0) weidis=0.1D0
1478 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1479 write (iout,'(a,i5)') 'nexcl: ',nexcl
1480 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1483 c------------------------------------------------------------------------------
1484 subroutine setup_var
1485 implicit real*8 (a-h,o-z)
1486 include 'DIMENSIONS'
1487 include 'COMMON.IOUNITS'
1488 include 'COMMON.GEO'
1489 include 'COMMON.VAR'
1490 include 'COMMON.INTERACT'
1491 include 'COMMON.LOCAL'
1492 include 'COMMON.NAMES'
1493 include 'COMMON.CHAIN'
1494 include 'COMMON.FFIELD'
1495 include 'COMMON.SBRIDGE'
1496 include 'COMMON.HEADER'
1497 include 'COMMON.CONTROL'
1498 include 'COMMON.DBASE'
1499 include 'COMMON.THREAD'
1500 include 'COMMON.TIME1'
1501 C Set up variable list.
1507 if (itype(i).ne.10) then
1509 ialph(i,1)=nvar+nside
1513 if (indphi.gt.0) then
1515 else if (indback.gt.0) then
1520 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1523 c----------------------------------------------------------------------------
1524 subroutine gen_dist_constr
1525 C Generate CA distance constraints.
1526 implicit real*8 (a-h,o-z)
1527 include 'DIMENSIONS'
1528 include 'COMMON.IOUNITS'
1529 include 'COMMON.GEO'
1530 include 'COMMON.VAR'
1531 include 'COMMON.INTERACT'
1532 include 'COMMON.LOCAL'
1533 include 'COMMON.NAMES'
1534 include 'COMMON.CHAIN'
1535 include 'COMMON.FFIELD'
1536 include 'COMMON.SBRIDGE'
1537 include 'COMMON.HEADER'
1538 include 'COMMON.CONTROL'
1539 include 'COMMON.DBASE'
1540 include 'COMMON.THREAD'
1541 include 'COMMON.TIME1'
1542 dimension itype_pdb(maxres)
1543 common /pizda/ itype_pdb
1545 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1546 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1547 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1549 do i=nstart_sup,nstart_sup+nsup-1
1550 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1551 cd & ' seq_pdb', restyp(itype_pdb(i))
1552 do j=i+2,nstart_sup+nsup-1
1554 ihpb(nhpb)=i+nstart_seq-nstart_sup
1555 jhpb(nhpb)=j+nstart_seq-nstart_sup
1557 dhpb(nhpb)=dist(i,j)
1560 cd write (iout,'(a)') 'Distance constraints:'
1565 cd if (ii.gt.nres) then
1570 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1571 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1572 cd & dhpb(i),forcon(i)
1576 c----------------------------------------------------------------------------
1578 implicit real*8 (a-h,o-z)
1579 include 'DIMENSIONS'
1580 include 'COMMON.MAP'
1581 include 'COMMON.IOUNITS'
1582 character*3 angid(4) /'THE','PHI','ALP','OME'/
1583 character*80 mapcard,ucase
1585 read (inp,'(a)') mapcard
1586 mapcard=ucase(mapcard)
1587 if (index(mapcard,'PHI').gt.0) then
1589 else if (index(mapcard,'THE').gt.0) then
1591 else if (index(mapcard,'ALP').gt.0) then
1593 else if (index(mapcard,'OME').gt.0) then
1596 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1597 stop 'Error - illegal variable spec in MAP card.'
1599 call readi (mapcard,'RES1',res1(imap),0)
1600 call readi (mapcard,'RES2',res2(imap),0)
1601 if (res1(imap).eq.0) then
1602 res1(imap)=res2(imap)
1603 else if (res2(imap).eq.0) then
1604 res2(imap)=res1(imap)
1606 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1608 & 'Error - illegal definition of variable group in MAP.'
1609 stop 'Error - illegal definition of variable group in MAP.'
1611 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1612 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1613 call readi(mapcard,'NSTEP',nstep(imap),0)
1614 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1616 & 'Illegal boundary and/or step size specification in MAP.'
1617 stop 'Illegal boundary and/or step size specification in MAP.'
1622 c----------------------------------------------------------------------------
1623 csa subroutine csaread
1624 csa implicit real*8 (a-h,o-z)
1625 csa include 'DIMENSIONS'
1626 csa include 'COMMON.IOUNITS'
1627 csa include 'COMMON.GEO'
1628 csa include 'COMMON.CSA'
1629 csa include 'COMMON.BANK'
1630 csa include 'COMMON.CONTROL'
1631 csa character*80 ucase
1632 csa character*620 mcmcard
1633 csa call card_concat(mcmcard)
1635 csa call readi(mcmcard,'NCONF',nconf,50)
1636 csa call readi(mcmcard,'NADD',nadd,0)
1637 csa call readi(mcmcard,'JSTART',jstart,1)
1638 csa call readi(mcmcard,'JEND',jend,1)
1639 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1640 csa call readi(mcmcard,'N0',n0,1)
1641 csa call readi(mcmcard,'N1',n1,6)
1642 csa call readi(mcmcard,'N2',n2,4)
1643 csa call readi(mcmcard,'N3',n3,0)
1644 csa call readi(mcmcard,'N4',n4,0)
1645 csa call readi(mcmcard,'N5',n5,0)
1646 csa call readi(mcmcard,'N6',n6,10)
1647 csa call readi(mcmcard,'N7',n7,0)
1648 csa call readi(mcmcard,'N8',n8,0)
1649 csa call readi(mcmcard,'N9',n9,0)
1650 csa call readi(mcmcard,'N14',n14,0)
1651 csa call readi(mcmcard,'N15',n15,0)
1652 csa call readi(mcmcard,'N16',n16,0)
1653 csa call readi(mcmcard,'N17',n17,0)
1654 csa call readi(mcmcard,'N18',n18,0)
1656 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1658 csa call readi(mcmcard,'NDIFF',ndiff,2)
1659 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1660 csa call readi(mcmcard,'IS1',is1,1)
1661 csa call readi(mcmcard,'IS2',is2,8)
1662 csa call readi(mcmcard,'NRAN0',nran0,4)
1663 csa call readi(mcmcard,'NRAN1',nran1,2)
1664 csa call readi(mcmcard,'IRR',irr,1)
1665 csa call readi(mcmcard,'NSEED',nseed,20)
1666 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1667 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1668 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1669 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1670 csa call readi(mcmcard,'ICMAX',icmax,3)
1671 csa call readi(mcmcard,'IRESTART',irestart,0)
1672 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1675 csa call reada(mcmcard,'DELE',dele,20.0d0)
1676 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1677 csa call readi(mcmcard,'IREF',iref,0)
1678 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1679 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1680 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1681 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1682 csa write (iout,*) "NCONF_IN",nconf_in
1685 c----------------------------------------------------------------------------
1686 cfmc subroutine mcmfread
1687 cfmc implicit real*8 (a-h,o-z)
1688 cfmc include 'DIMENSIONS'
1689 cfmc include 'COMMON.MCMF'
1690 cfmc include 'COMMON.IOUNITS'
1691 cfmc include 'COMMON.GEO'
1692 cfmc character*80 ucase
1693 cfmc character*620 mcmcard
1694 cfmc call card_concat(mcmcard)
1696 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1697 cfmc write(iout,*)'MAXRANT=',maxrant
1698 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1699 cfmc write(iout,*)'MAXFAM=',maxfam
1700 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1701 cfmc write(iout,*)'NNET1=',nnet1
1702 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1703 cfmc write(iout,*)'NNET2=',nnet2
1704 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1705 cfmc write(iout,*)'NNET3=',nnet3
1706 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1707 cfmc write(iout,*)'ILASTT=',ilastt
1708 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1709 cfmc write(iout,*)'MAXSTR=',maxstr
1710 cfmc maxstr_f=maxstr/maxfam
1711 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1712 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1713 cfmc write(iout,*)'NMCMF=',nmcmf
1714 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1715 cfmc write(iout,*)'IFOCUS=',ifocus
1716 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1717 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1718 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1719 cfmc write(iout,*)'INTPRT=',intprt
1720 cfmc call readi(mcmcard,'IPRT',iprt,100)
1721 cfmc write(iout,*)'IPRT=',iprt
1722 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1723 cfmc write(iout,*)'IMAXTR=',imaxtr
1724 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1725 cfmc write(iout,*)'MAXEVEN=',maxeven
1726 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1727 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1728 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1729 cfmc write(iout,*)'INIMIN=',inimin
1730 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1731 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1732 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1733 cfmc write(iout,*)'NTHREAD=',nthread
1734 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1735 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1736 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1737 cfmc write(iout,*)'MAXPERT=',maxpert
1738 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1739 cfmc write(iout,*)'IRMSD=',irmsd
1740 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1741 cfmc write(iout,*)'DENEMIN=',denemin
1742 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1743 cfmc write(iout,*)'RCUT1S=',rcut1s
1744 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1745 cfmc write(iout,*)'RCUT1E=',rcut1e
1746 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1747 cfmc write(iout,*)'RCUT2S=',rcut2s
1748 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1749 cfmc write(iout,*)'RCUT2E=',rcut2e
1750 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1751 cfmc write(iout,*)'DPERT1=',d_pert1
1752 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1753 cfmc write(iout,*)'DPERT1A=',d_pert1a
1754 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1755 cfmc write(iout,*)'DPERT2=',d_pert2
1756 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1757 cfmc write(iout,*)'DPERT2A=',d_pert2a
1758 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1759 cfmc write(iout,*)'DPERT2B=',d_pert2b
1760 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1761 cfmc write(iout,*)'DPERT2C=',d_pert2c
1762 cfmc d_pert1=deg2rad*d_pert1
1763 cfmc d_pert1a=deg2rad*d_pert1a
1764 cfmc d_pert2=deg2rad*d_pert2
1765 cfmc d_pert2a=deg2rad*d_pert2a
1766 cfmc d_pert2b=deg2rad*d_pert2b
1767 cfmc d_pert2c=deg2rad*d_pert2c
1768 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1769 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1770 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1771 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1772 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1773 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1774 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1775 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1776 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1777 cfmc write(iout,*)'RCUTINI=',rcutini
1778 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1779 cfmc write(iout,*)'GRAT=',grat
1780 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1781 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1785 c----------------------------------------------------------------------------
1787 implicit real*8 (a-h,o-z)
1788 include 'DIMENSIONS'
1789 include 'COMMON.MCM'
1790 include 'COMMON.MCE'
1791 include 'COMMON.IOUNITS'
1793 character*320 mcmcard
1794 call card_concat(mcmcard)
1795 call readi(mcmcard,'MAXACC',maxacc,100)
1796 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1797 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1798 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1799 call readi(mcmcard,'MAXREPM',maxrepm,200)
1800 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1801 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1802 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1803 call reada(mcmcard,'E_UP',e_up,5.0D0)
1804 call reada(mcmcard,'DELTE',delte,0.1D0)
1805 call readi(mcmcard,'NSWEEP',nsweep,5)
1806 call readi(mcmcard,'NSTEPH',nsteph,0)
1807 call readi(mcmcard,'NSTEPC',nstepc,0)
1808 call reada(mcmcard,'TMIN',tmin,298.0D0)
1809 call reada(mcmcard,'TMAX',tmax,298.0D0)
1810 call readi(mcmcard,'NWINDOW',nwindow,0)
1811 call readi(mcmcard,'PRINT_MC',print_mc,0)
1812 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1813 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1814 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1815 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1816 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1817 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1818 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1819 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1820 if (nwindow.gt.0) then
1821 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1823 winlen(i)=winend(i)-winstart(i)+1
1826 if (tmax.lt.tmin) tmax=tmin
1827 if (tmax.eq.tmin) then
1831 if (nstepc.gt.0 .and. nsteph.gt.0) then
1832 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1833 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1835 C Probabilities of different move types
1836 sumpro_type(0)=0.0D0
1837 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1838 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1839 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1840 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1841 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1842 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1843 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1845 print *,'i',i,' sumprotype',sumpro_type(i)
1846 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1847 print *,'i',i,' sumprotype',sumpro_type(i)
1851 c----------------------------------------------------------------------------
1852 subroutine read_minim
1853 implicit real*8 (a-h,o-z)
1854 include 'DIMENSIONS'
1855 include 'COMMON.MINIM'
1856 include 'COMMON.IOUNITS'
1858 character*320 minimcard
1859 call card_concat(minimcard)
1860 call readi(minimcard,'MAXMIN',maxmin,2000)
1861 call readi(minimcard,'MAXFUN',maxfun,5000)
1862 call readi(minimcard,'MINMIN',minmin,maxmin)
1863 call readi(minimcard,'MINFUN',minfun,maxmin)
1864 call reada(minimcard,'TOLF',tolf,1.0D-2)
1865 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1866 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1867 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1868 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1869 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1870 & 'Options in energy minimization:'
1871 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1872 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1873 & 'MinMin:',MinMin,' MinFun:',MinFun,
1874 & ' TolF:',TolF,' RTolF:',RTolF
1877 c----------------------------------------------------------------------------
1878 subroutine read_angles(kanal,*)
1879 implicit real*8 (a-h,o-z)
1880 include 'DIMENSIONS'
1881 include 'COMMON.GEO'
1882 include 'COMMON.VAR'
1883 include 'COMMON.CHAIN'
1884 include 'COMMON.IOUNITS'
1885 include 'COMMON.CONTROL'
1886 c Read angles from input
1888 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1889 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1890 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1891 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1894 c 9/7/01 avoid 180 deg valence angle
1895 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1897 theta(i)=deg2rad*theta(i)
1898 phi(i)=deg2rad*phi(i)
1899 alph(i)=deg2rad*alph(i)
1900 omeg(i)=deg2rad*omeg(i)
1905 c----------------------------------------------------------------------------
1906 subroutine reada(rekord,lancuch,wartosc,default)
1908 character*(*) rekord,lancuch
1909 double precision wartosc,default
1912 iread=index(rekord,lancuch)
1913 if (iread.eq.0) then
1917 iread=iread+ilen(lancuch)+1
1918 read (rekord(iread:),*,err=10,end=10) wartosc
1923 c----------------------------------------------------------------------------
1924 subroutine readi(rekord,lancuch,wartosc,default)
1926 character*(*) rekord,lancuch
1927 integer wartosc,default
1930 iread=index(rekord,lancuch)
1931 if (iread.eq.0) then
1935 iread=iread+ilen(lancuch)+1
1936 read (rekord(iread:),*,err=10,end=10) wartosc
1941 c----------------------------------------------------------------------------
1942 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1945 integer tablica(dim),default
1946 character*(*) rekord,lancuch
1953 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1954 if (iread.eq.0) return
1955 iread=iread+ilen(lancuch)+1
1956 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1959 c----------------------------------------------------------------------------
1960 subroutine multreada(rekord,lancuch,tablica,dim,default)
1963 double precision tablica(dim),default
1964 character*(*) rekord,lancuch
1971 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1972 if (iread.eq.0) return
1973 iread=iread+ilen(lancuch)+1
1974 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1977 c----------------------------------------------------------------------------
1978 subroutine openunits
1979 implicit real*8 (a-h,o-z)
1980 include 'DIMENSIONS'
1983 character*16 form,nodename
1986 include 'COMMON.SETUP'
1987 include 'COMMON.IOUNITS'
1989 include 'COMMON.CONTROL'
1990 integer lenpre,lenpot,ilen,lentmp
1992 character*3 out1file_text,ucase
1995 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1996 call getenv_loc("PREFIX",prefix)
1998 call getenv_loc("POT",pot)
1999 call getenv_loc("DIRTMP",tmpdir)
2000 call getenv_loc("CURDIR",curdir)
2001 call getenv_loc("OUT1FILE",out1file_text)
2002 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2003 out1file_text=ucase(out1file_text)
2004 if (out1file_text(1:1).eq."Y") then
2007 out1file=fg_rank.gt.0
2012 if (lentmp.gt.0) then
2013 write (*,'(80(1h!))')
2014 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2015 write (*,'(80(1h!))')
2016 write (*,*)"All output files will be on node /tmp directory."
2018 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2019 if (me.eq.king) then
2020 write (*,*) "The master node is ",nodename
2021 else if (fg_rank.eq.0) then
2022 write (*,*) "I am the CG slave node ",nodename
2024 write (*,*) "I am the FG slave node ",nodename
2027 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2028 lenpre = lentmp+lenpre+1
2030 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2031 C Get the names and open the input files
2032 #if defined(WINIFL) || defined(WINPGI)
2033 open(1,file=pref_orig(:ilen(pref_orig))//
2034 & '.inp',status='old',readonly,shared)
2035 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2036 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2037 C Get parameter filenames and open the parameter files.
2038 call getenv_loc('BONDPAR',bondname)
2039 open (ibond,file=bondname,status='old',readonly,shared)
2040 call getenv_loc('THETPAR',thetname)
2041 open (ithep,file=thetname,status='old',readonly,shared)
2043 call getenv_loc('THETPARPDB',thetname_pdb)
2044 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2046 call getenv_loc('ROTPAR',rotname)
2047 open (irotam,file=rotname,status='old',readonly,shared)
2049 call getenv_loc('ROTPARPDB',rotname_pdb)
2050 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2052 call getenv_loc('TORPAR',torname)
2053 open (itorp,file=torname,status='old',readonly,shared)
2054 call getenv_loc('TORDPAR',tordname)
2055 open (itordp,file=tordname,status='old',readonly,shared)
2056 call getenv_loc('FOURIER',fouriername)
2057 open (ifourier,file=fouriername,status='old',readonly,shared)
2058 call getenv_loc('ELEPAR',elename)
2059 open (ielep,file=elename,status='old',readonly,shared)
2060 call getenv_loc('SIDEPAR',sidename)
2061 open (isidep,file=sidename,status='old',readonly,shared)
2062 #elif (defined CRAY) || (defined AIX)
2063 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2065 c print *,"Processor",myrank," opened file 1"
2066 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2067 c print *,"Processor",myrank," opened file 9"
2068 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2069 C Get parameter filenames and open the parameter files.
2070 call getenv_loc('BONDPAR',bondname)
2071 open (ibond,file=bondname,status='old',action='read')
2072 c print *,"Processor",myrank," opened file IBOND"
2073 call getenv_loc('THETPAR',thetname)
2074 open (ithep,file=thetname,status='old',action='read')
2075 c print *,"Processor",myrank," opened file ITHEP"
2077 call getenv_loc('THETPARPDB',thetname_pdb)
2078 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2080 call getenv_loc('ROTPAR',rotname)
2081 open (irotam,file=rotname,status='old',action='read')
2082 c print *,"Processor",myrank," opened file IROTAM"
2084 call getenv_loc('ROTPARPDB',rotname_pdb)
2085 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2087 call getenv_loc('TORPAR',torname)
2088 open (itorp,file=torname,status='old',action='read')
2089 c print *,"Processor",myrank," opened file ITORP"
2090 call getenv_loc('TORDPAR',tordname)
2091 open (itordp,file=tordname,status='old',action='read')
2092 c print *,"Processor",myrank," opened file ITORDP"
2093 call getenv_loc('SCCORPAR',sccorname)
2094 open (isccor,file=sccorname,status='old',action='read')
2095 c print *,"Processor",myrank," opened file ISCCOR"
2096 call getenv_loc('FOURIER',fouriername)
2097 open (ifourier,file=fouriername,status='old',action='read')
2098 c print *,"Processor",myrank," opened file IFOURIER"
2099 call getenv_loc('ELEPAR',elename)
2100 open (ielep,file=elename,status='old',action='read')
2101 c print *,"Processor",myrank," opened file IELEP"
2102 call getenv_loc('SIDEPAR',sidename)
2103 open (isidep,file=sidename,status='old',action='read')
2104 c print *,"Processor",myrank," opened file ISIDEP"
2105 c print *,"Processor",myrank," opened parameter files"
2107 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2108 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2109 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2110 C Get parameter filenames and open the parameter files.
2111 call getenv_loc('BONDPAR',bondname)
2112 open (ibond,file=bondname,status='old')
2113 call getenv_loc('THETPAR',thetname)
2114 open (ithep,file=thetname,status='old')
2116 call getenv_loc('THETPARPDB',thetname_pdb)
2117 open (ithep_pdb,file=thetname_pdb,status='old')
2119 call getenv_loc('ROTPAR',rotname)
2120 open (irotam,file=rotname,status='old')
2122 call getenv_loc('ROTPARPDB',rotname_pdb)
2123 open (irotam_pdb,file=rotname_pdb,status='old')
2125 call getenv_loc('TORPAR',torname)
2126 open (itorp,file=torname,status='old')
2127 call getenv_loc('TORDPAR',tordname)
2128 open (itordp,file=tordname,status='old')
2129 call getenv_loc('SCCORPAR',sccorname)
2130 open (isccor,file=sccorname,status='old')
2131 call getenv_loc('FOURIER',fouriername)
2132 open (ifourier,file=fouriername,status='old')
2133 call getenv_loc('ELEPAR',elename)
2134 open (ielep,file=elename,status='old')
2135 call getenv_loc('SIDEPAR',sidename)
2136 open (isidep,file=sidename,status='old')
2138 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2140 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2141 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2142 C Get parameter filenames and open the parameter files.
2143 call getenv_loc('BONDPAR',bondname)
2144 open (ibond,file=bondname,status='old',action='read')
2145 call getenv_loc('THETPAR',thetname)
2146 open (ithep,file=thetname,status='old',action='read')
2148 call getenv_loc('THETPARPDB',thetname_pdb)
2149 print *,"thetname_pdb ",thetname_pdb
2150 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2151 print *,ithep_pdb," opened"
2153 call getenv_loc('ROTPAR',rotname)
2154 open (irotam,file=rotname,status='old',action='read')
2156 call getenv_loc('ROTPARPDB',rotname_pdb)
2157 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2159 call getenv_loc('TORPAR',torname)
2160 open (itorp,file=torname,status='old',action='read')
2161 call getenv_loc('TORDPAR',tordname)
2162 open (itordp,file=tordname,status='old',action='read')
2163 call getenv_loc('SCCORPAR',sccorname)
2164 open (isccor,file=sccorname,status='old',action='read')
2165 call getenv_loc('FOURIER',fouriername)
2166 open (ifourier,file=fouriername,status='old',action='read')
2167 call getenv_loc('ELEPAR',elename)
2168 open (ielep,file=elename,status='old',action='read')
2169 call getenv_loc('SIDEPAR',sidename)
2170 open (isidep,file=sidename,status='old',action='read')
2174 C 8/9/01 In the newest version SCp interaction constants are read from a file
2175 C Use -DOLDSCP to use hard-coded constants instead.
2177 call getenv_loc('SCPPAR',scpname)
2178 #if defined(WINIFL) || defined(WINPGI)
2179 open (iscpp,file=scpname,status='old',readonly,shared)
2180 #elif (defined CRAY) || (defined AIX)
2181 open (iscpp,file=scpname,status='old',action='read')
2183 open (iscpp,file=scpname,status='old')
2185 open (iscpp,file=scpname,status='old',action='read')
2188 call getenv_loc('PATTERN',patname)
2189 #if defined(WINIFL) || defined(WINPGI)
2190 open (icbase,file=patname,status='old',readonly,shared)
2191 #elif (defined CRAY) || (defined AIX)
2192 open (icbase,file=patname,status='old',action='read')
2194 open (icbase,file=patname,status='old')
2196 open (icbase,file=patname,status='old',action='read')
2199 C Open output file only for CG processes
2200 c print *,"Processor",myrank," fg_rank",fg_rank
2201 if (fg_rank.eq.0) then
2203 if (nodes.eq.1) then
2206 npos = dlog10(dfloat(nodes-1))+1
2208 if (npos.lt.3) npos=3
2209 write (liczba,'(i1)') npos
2210 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2212 write (liczba,form) me
2213 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2214 & liczba(:ilen(liczba))
2215 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2217 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2219 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2220 & liczba(:ilen(liczba))//'.mol2'
2221 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2222 & liczba(:ilen(liczba))//'.stat'
2224 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2225 & //liczba(:ilen(liczba))//'.stat')
2226 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2229 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2230 & liczba(:ilen(liczba))//'.const'
2235 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2236 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2237 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2238 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2239 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2241 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2243 rest2name=prefix(:ilen(prefix))//'.rst'
2245 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2248 #if defined(AIX) || defined(PGI)
2249 if (me.eq.king .or. .not. out1file)
2250 & open(iout,file=outname,status='unknown')
2253 if (fg_rank.gt.0) then
2254 write (liczba,'(i3.3)') myrank/nfgtasks
2255 write (ll,'(bz,i3.3)') fg_rank
2256 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2262 open(igeom,file=intname,status='unknown',position='append')
2263 open(ipdb,file=pdbname,status='unknown')
2264 open(imol2,file=mol2name,status='unknown')
2265 open(istat,file=statname,status='unknown',position='append')
2267 c1out open(iout,file=outname,status='unknown')
2270 if (me.eq.king .or. .not.out1file)
2271 & open(iout,file=outname,status='unknown')
2274 if (fg_rank.gt.0) then
2275 print "Processor",fg_rank," opening output file"
2276 write (liczba,'(i3.3)') myrank/nfgtasks
2277 write (ll,'(bz,i3.3)') fg_rank
2278 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2284 open(igeom,file=intname,status='unknown',access='append')
2285 open(ipdb,file=pdbname,status='unknown')
2286 open(imol2,file=mol2name,status='unknown')
2287 open(istat,file=statname,status='unknown',access='append')
2289 c1out open(iout,file=outname,status='unknown')
2292 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2293 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2294 csa csa_history=prefix(:lenpre)//'.CSA.history'
2295 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2296 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2297 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2298 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2299 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2300 csa csa_int=prefix(:lenpre)//'.int'
2301 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2302 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2303 csa csa_in=prefix(:lenpre)//'.CSA.in'
2304 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2307 write (iout,'(80(1h-))')
2308 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2309 write (iout,'(80(1h-))')
2310 write (iout,*) "Input file : ",
2311 & pref_orig(:ilen(pref_orig))//'.inp'
2312 write (iout,*) "Output file : ",
2313 & outname(:ilen(outname))
2315 write (iout,*) "Sidechain potential file : ",
2316 & sidename(:ilen(sidename))
2318 write (iout,*) "SCp potential file : ",
2319 & scpname(:ilen(scpname))
2321 write (iout,*) "Electrostatic potential file : ",
2322 & elename(:ilen(elename))
2323 write (iout,*) "Cumulant coefficient file : ",
2324 & fouriername(:ilen(fouriername))
2325 write (iout,*) "Torsional parameter file : ",
2326 & torname(:ilen(torname))
2327 write (iout,*) "Double torsional parameter file : ",
2328 & tordname(:ilen(tordname))
2329 write (iout,*) "SCCOR parameter file : ",
2330 & sccorname(:ilen(sccorname))
2331 write (iout,*) "Bond & inertia constant file : ",
2332 & bondname(:ilen(bondname))
2333 write (iout,*) "Bending parameter file : ",
2334 & thetname(:ilen(thetname))
2335 write (iout,*) "Rotamer parameter file : ",
2336 & rotname(:ilen(rotname))
2337 write (iout,*) "Threading database : ",
2338 & patname(:ilen(patname))
2340 &write (iout,*)" DIRTMP : ",
2342 write (iout,'(80(1h-))')
2346 c----------------------------------------------------------------------------
2347 subroutine card_concat(card)
2348 implicit real*8 (a-h,o-z)
2349 include 'DIMENSIONS'
2350 include 'COMMON.IOUNITS'
2352 character*80 karta,ucase
2354 read (inp,'(a)') karta
2357 do while (karta(80:80).eq.'&')
2358 card=card(:ilen(card)+1)//karta(:79)
2359 read (inp,'(a)') karta
2362 card=card(:ilen(card)+1)//karta
2365 c----------------------------------------------------------------------------------
2367 implicit real*8 (a-h,o-z)
2368 include 'DIMENSIONS'
2369 include 'COMMON.CHAIN'
2370 include 'COMMON.IOUNITS'
2372 open(irest2,file=rest2name,status='unknown')
2373 read(irest2,*) totT,EK,potE,totE,t_bath
2375 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2378 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2381 read (irest2,*) iset
2386 c---------------------------------------------------------------------------------
2387 subroutine read_fragments
2388 implicit real*8 (a-h,o-z)
2389 include 'DIMENSIONS'
2393 include 'COMMON.SETUP'
2394 include 'COMMON.CHAIN'
2395 include 'COMMON.IOUNITS'
2397 include 'COMMON.CONTROL'
2398 read(inp,*) nset,nfrag,npair,nfrag_back
2399 if(me.eq.king.or..not.out1file)
2400 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2401 & " nfrag_back",nfrag_back
2403 read(inp,*) mset(iset)
2405 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2407 if(me.eq.king.or..not.out1file)
2408 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2409 & ifrag(2,i,iset), qinfrag(i,iset)
2412 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2414 if(me.eq.king.or..not.out1file)
2415 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2416 & ipair(2,i,iset), qinpair(i,iset)
2419 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2420 & wfrag_back(3,i,iset),
2421 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2422 if(me.eq.king.or..not.out1file)
2423 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2424 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2429 c-------------------------------------------------------------------------------
2430 subroutine read_dist_constr
2431 implicit real*8 (a-h,o-z)
2432 include 'DIMENSIONS'
2436 include 'COMMON.SETUP'
2437 include 'COMMON.CONTROL'
2438 include 'COMMON.CHAIN'
2439 include 'COMMON.IOUNITS'
2440 include 'COMMON.SBRIDGE'
2441 integer ifrag_(2,100),ipair_(2,100)
2442 double precision wfrag_(100),wpair_(100)
2443 character*500 controlcard
2444 c write (iout,*) "Calling read_dist_constr"
2445 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2447 call card_concat(controlcard)
2448 call readi(controlcard,"NFRAG",nfrag_,0)
2449 call readi(controlcard,"NPAIR",npair_,0)
2450 call readi(controlcard,"NDIST",ndist_,0)
2451 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2452 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2453 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2454 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2455 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2456 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2457 c write (iout,*) "IFRAG"
2459 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2461 c write (iout,*) "IPAIR"
2463 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2465 if (.not.refstr .and. nfrag.gt.0) then
2467 & "ERROR: no reference structure to compute distance restraints"
2469 & "Restraints must be specified explicitly (NDIST=number)"
2472 if (nfrag.lt.2 .and. npair.gt.0) then
2473 write (iout,*) "ERROR: Less than 2 fragments specified",
2474 & " but distance restraints between pairs requested"
2479 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2480 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2481 & ifrag_(2,i)=nstart_sup+nsup-1
2482 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2484 if (wfrag_(i).gt.0.0d0) then
2485 do j=ifrag_(1,i),ifrag_(2,i)-1
2486 do k=j+1,ifrag_(2,i)
2487 write (iout,*) "j",j," k",k
2489 if (constr_dist.eq.1) then
2494 forcon(nhpb)=wfrag_(i)
2495 else if (constr_dist.eq.2) then
2496 if (ddjk.le.dist_cut) then
2501 forcon(nhpb)=wfrag_(i)
2508 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2511 if (.not.out1file .or. me.eq.king)
2512 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2513 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2515 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2516 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2523 if (wpair_(i).gt.0.0d0) then
2531 do j=ifrag_(1,ii),ifrag_(2,ii)
2532 do k=ifrag_(1,jj),ifrag_(2,jj)
2536 forcon(nhpb)=wpair_(i)
2537 dhpb(nhpb)=dist(j,k)
2539 if (.not.out1file .or. me.eq.king)
2540 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2541 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2543 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2544 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2551 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2552 & ibecarb(i),forcon(nhpb+1)
2553 if (forcon(nhpb+1).gt.0.0d0) then
2555 if (ibecarb(i).gt.0) then
2556 ihpb(i)=ihpb(i)+nres
2557 jhpb(i)=jhpb(i)+nres
2559 if (dhpb(nhpb).eq.0.0d0)
2560 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2564 if (.not.out1file .or. me.eq.king) then
2567 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2568 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2576 c-------------------------------------------------------------------------------
2578 subroutine flush(iu)
2583 subroutine flush(iu)
2588 c------------------------------------------------------------------------------
2589 subroutine copy_to_tmp(source)
2590 include "DIMENSIONS"
2591 include "COMMON.IOUNITS"
2592 character*(*) source
2593 character* 256 tmpfile
2597 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2598 inquire(file=tmpfile,exist=ex)
2600 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2601 & " to temporary directory..."
2602 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2603 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2607 c------------------------------------------------------------------------------
2608 subroutine move_from_tmp(source)
2609 include "DIMENSIONS"
2610 include "COMMON.IOUNITS"
2611 character*(*) source
2614 write (*,*) "Moving ",source(:ilen(source)),
2615 & " from temporary directory to working directory"
2616 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2617 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2620 c------------------------------------------------------------------------------
2621 subroutine random_init(seed)
2623 C Initialize random number generator
2625 implicit real*8 (a-h,o-z)
2626 include 'DIMENSIONS'
2632 logical OKRandom, prng_restart
2634 integer iseed_array(4)
2636 include 'COMMON.IOUNITS'
2637 include 'COMMON.TIME1'
2638 include 'COMMON.THREAD'
2639 include 'COMMON.SBRIDGE'
2640 include 'COMMON.CONTROL'
2641 include 'COMMON.MCM'
2642 include 'COMMON.MAP'
2643 include 'COMMON.HEADER'
2644 csa include 'COMMON.CSA'
2645 include 'COMMON.CHAIN'
2646 include 'COMMON.MUCA'
2648 include 'COMMON.FFIELD'
2649 include 'COMMON.SETUP'
2650 iseed=-dint(dabs(seed))
2651 if (iseed.eq.0) then
2652 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2653 & 'Random seed undefined. The program will stop.'
2654 write (*,'(/80(1h*)/20x,a/80(1h*))')
2655 & 'Random seed undefined. The program will stop.'
2657 call mpi_finalize(mpi_comm_world,ierr)
2659 stop 'Bad random seed.'
2662 if (fg_rank.eq.0) then
2666 if(me.eq.king .or. .not. out1file)
2667 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2668 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2669 OKRandom = prng_restart(me,iseedi8)
2672 tmp=65536.0d0**(4-i)
2673 iseed_array(i) = dint(seed/tmp)
2674 seed=seed-iseed_array(i)*tmp
2676 if(me.eq.king .or. .not. out1file)
2677 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2678 & (iseed_array(i),i=1,4)
2679 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2680 & (iseed_array(i),i=1,4)
2681 OKRandom = prng_restart(me,iseed_array)
2684 r1=ran_number(0.0D0,1.0D0)
2685 if(me.eq.king .or. .not. out1file)
2686 & write (iout,*) 'ran_num',r1
2687 if (r1.lt.0.0d0) OKRandom=.false.
2689 if (.not.OKRandom) then
2690 write (iout,*) 'PRNG IS NOT WORKING!!!'
2691 print *,'PRNG IS NOT WORKING!!!'
2694 call mpi_abort(mpi_comm_world,error_msg,ierr)
2697 write (iout,*) 'too many processors for parallel prng'
2698 write (*,*) 'too many processors for parallel prng'
2706 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)