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 c 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 if (itype(2).eq.10) then
961 & "Glycine is the first full residue, initial dummy deleted"
967 if (itype(nres).eq.10) then
969 & "Glycine is the last full residue, terminal dummy deleted"
973 C Assign initial virtual bond lengths
979 vbld(i+nres)=dsc(itype(i))
980 vbld_inv(i+nres)=dsc_inv(itype(i))
981 c write (iout,*) "i",i," itype",itype(i),
982 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
986 c print '(20i4)',(itype(i),i=1,nres)
989 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
991 if (itype(i).eq.21) then
995 else if (itype(i+1).ne.20) then
997 else if (itype(i).ne.20) then
1004 if(me.eq.king.or..not.out1file)then
1005 write (iout,*) "ITEL"
1007 write (iout,*) i,itype(i),itel(i)
1009 print *,'Call Read_Bridge.'
1012 C 8/13/98 Set limits to generating the dihedral angles
1017 read (inp,*) ndih_constr
1018 if (ndih_constr.gt.0) then
1020 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1021 if(me.eq.king.or..not.out1file)then
1023 & 'There are',ndih_constr,' constraints on phi angles.'
1025 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1029 phi0(i)=deg2rad*phi0(i)
1030 drange(i)=deg2rad*drange(i)
1032 if(me.eq.king.or..not.out1file)
1033 & write (iout,*) 'FTORS',ftors
1036 phibound(1,ii) = phi0(i)-drange(i)
1037 phibound(2,ii) = phi0(i)+drange(i)
1042 if (me.eq.king) then
1044 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1046 write (iout,'(a3,i5,2f10.1)')
1047 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1053 cd print *,'NNT=',NNT,' NCT=',NCT
1054 if (itype(1).eq.21) nnt=2
1055 if (itype(nres).eq.21) nct=nct-1
1057 if(me.eq.king.or..not.out1file)
1058 & write (iout,'(a,i3)') 'nsup=',nsup
1060 if (nsup.le.(nct-nnt+1)) then
1061 do i=0,nct-nnt+1-nsup
1062 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1068 & 'Error - sequences to be superposed do not match.'
1071 do i=0,nsup-(nct-nnt+1)
1072 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1074 nstart_sup=nstart_sup+i
1080 & 'Error - sequences to be superposed do not match.'
1083 if (nsup.eq.0) nsup=nct-nnt
1084 if (nstart_sup.eq.0) nstart_sup=nnt
1085 if (nstart_seq.eq.0) nstart_seq=nnt
1086 if(me.eq.king.or..not.out1file)
1087 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1088 & ' nstart_seq=',nstart_seq
1090 c--- Zscore rms -------
1091 if (nz_start.eq.0) nz_start=nnt
1092 if (nz_end.eq.0 .and. nsup.gt.0) then
1094 else if (nz_end.eq.0) then
1097 if(me.eq.king.or..not.out1file)then
1098 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1099 write (iout,*) 'IZ_SC=',iz_sc
1101 c----------------------
1104 if (.not.pdbref) then
1105 call read_angles(inp,*38)
1107 38 write (iout,'(a)') 'Error reading reference structure.'
1109 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1110 stop 'Error reading reference structure'
1114 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1123 call contact(.true.,ncont_ref,icont_ref,co)
1125 if(me.eq.king.or..not.out1file)
1126 & write (iout,*) 'Contact order:',co
1128 if(me.eq.king.or..not.out1file)
1129 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1132 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1134 if(me.eq.king.or..not.out1file)
1135 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1136 & icont_ref(1,i),' ',
1137 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1141 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1142 if (constr_dist.gt.0) then
1143 call read_dist_constr
1145 if (nhpb.gt.0) call hpb_partition
1146 c write (iout,*) "After read_dist_constr nhpb",nhpb
1148 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1149 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1150 & modecalc.ne.10) then
1151 C If input structure hasn't been supplied from the PDB file read or generate
1153 if (iranconf.eq.0 .and. .not. extconf) then
1154 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1155 & write (iout,'(a)') 'Initial geometry will be read in.'
1157 read(inp,'(8f10.5)',end=36,err=36)
1158 & ((c(l,k),l=1,3),k=1,nres),
1159 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1160 call int_from_cart1(.false.)
1163 dc(j,i)=c(j,i+1)-c(j,i)
1164 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1168 if (itype(i).ne.10) then
1170 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1171 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1177 call read_angles(inp,*36)
1180 36 write (iout,'(a)') 'Error reading angle file.'
1182 call mpi_finalize( MPI_COMM_WORLD,IERR )
1184 stop 'Error reading angle file.'
1186 else if (extconf) then
1187 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1188 & write (iout,'(a)') 'Extended chain initial geometry.'
1190 theta(i)=90d0*deg2rad
1193 phi(i)=180d0*deg2rad
1196 alph(i)=110d0*deg2rad
1199 omeg(i)=-120d0*deg2rad
1202 if(me.eq.king.or..not.out1file)
1203 & write (iout,'(a)') 'Random-generated initial geometry.'
1207 if (me.eq.king .or. fg_rank.eq.0 .and. (
1208 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1212 call gen_rand_conf(itmp,*30)
1214 30 write (iout,*) 'Failed to generate random conformation',
1215 & ', itrial=',itrial
1216 write (*,*) 'Processor:',me,
1217 & ' Failed to generate random conformation',
1227 write (iout,'(a,i3,a)') 'Processor:',me,
1228 & ' error in generating random conformation.'
1229 write (*,'(a,i3,a)') 'Processor:',me,
1230 & ' error in generating random conformation.'
1233 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1240 elseif (modecalc.eq.4) then
1241 read (inp,'(a)') intinname
1242 open (intin,file=intinname,status='old',err=333)
1243 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1244 & write (iout,'(a)') 'intinname',intinname
1245 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1247 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1249 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1251 stop 'Error opening angle file.'
1255 C Generate distance constraints, if the PDB structure is to be regularized.
1256 if (nthread.gt.0) then
1257 call read_threadbase
1260 if (me.eq.king .or. .not. out1file)
1262 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1263 write (iout,'(/a,i3,a)')
1264 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1265 write (iout,'(20i4)') (iss(i),i=1,ns)
1267 write(iout,*)"Running with dynamic disulfide-bond formation"
1269 write (iout,'(/a/)') 'Pre-formed links are:'
1275 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1276 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1282 if (ns.gt.0.and.dyn_ss) then
1286 forcon(i-nss)=forcon(i)
1293 dyn_ss_mask(iss(i))=.true.
1296 if (i2ndstr.gt.0) call secstrp2dihc
1297 c call geom_to_var(nvar,x)
1298 c call etotal(energia(0))
1299 c call enerprint(energia(0))
1300 c call briefout(0,etot)
1302 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1303 cd write (iout,'(a)') 'Variable list:'
1304 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1306 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1307 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1308 & 'Processor',myrank,': end reading molecular data.'
1312 c--------------------------------------------------------------------------
1313 logical function seq_comp(itypea,itypeb,length)
1315 integer length,itypea(length),itypeb(length)
1318 if (itypea(i).ne.itypeb(i)) then
1326 c-----------------------------------------------------------------------------
1327 subroutine read_bridge
1328 C Read information about disulfide bridges.
1329 implicit real*8 (a-h,o-z)
1330 include 'DIMENSIONS'
1334 include 'COMMON.IOUNITS'
1335 include 'COMMON.GEO'
1336 include 'COMMON.VAR'
1337 include 'COMMON.INTERACT'
1338 include 'COMMON.LOCAL'
1339 include 'COMMON.NAMES'
1340 include 'COMMON.CHAIN'
1341 include 'COMMON.FFIELD'
1342 include 'COMMON.SBRIDGE'
1343 include 'COMMON.HEADER'
1344 include 'COMMON.CONTROL'
1345 include 'COMMON.DBASE'
1346 include 'COMMON.THREAD'
1347 include 'COMMON.TIME1'
1348 include 'COMMON.SETUP'
1349 C Read bridging residues.
1350 read (inp,*) ns,(iss(i),i=1,ns)
1352 if(me.eq.king.or..not.out1file)
1353 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1354 C Check whether the specified bridging residues are cystines.
1356 if (itype(iss(i)).ne.1) then
1357 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1358 & 'Do you REALLY think that the residue ',
1359 & restyp(itype(iss(i))),i,
1360 & ' can form a disulfide bridge?!!!'
1361 write (*,'(2a,i3,a)')
1362 & 'Do you REALLY think that the residue ',
1363 & restyp(itype(iss(i))),i,
1364 & ' can form a disulfide bridge?!!!'
1366 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1371 C Read preformed bridges.
1373 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1375 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1378 C Check if the residues involved in bridges are in the specified list of
1379 C bridging residues.
1382 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1383 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1384 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1385 & ' contains residues present in other pairs.'
1386 write (*,'(a,i3,a)') 'Disulfide pair',i,
1387 & ' contains residues present in other pairs.'
1389 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1395 if (ihpb(i).eq.iss(j)) goto 10
1397 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1400 if (jhpb(i).eq.iss(j)) goto 20
1402 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1408 ihpb(i)=ihpb(i)+nres
1409 jhpb(i)=jhpb(i)+nres
1415 c----------------------------------------------------------------------------
1416 subroutine read_x(kanal,*)
1417 implicit real*8 (a-h,o-z)
1418 include 'DIMENSIONS'
1419 include 'COMMON.GEO'
1420 include 'COMMON.VAR'
1421 include 'COMMON.CHAIN'
1422 include 'COMMON.IOUNITS'
1423 include 'COMMON.CONTROL'
1424 include 'COMMON.LOCAL'
1425 include 'COMMON.INTERACT'
1426 c Read coordinates from input
1428 read(kanal,'(8f10.5)',end=10,err=10)
1429 & ((c(l,k),l=1,3),k=1,nres),
1430 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1433 c(j,2*nres)=c(j,nres)
1435 call int_from_cart1(.false.)
1438 dc(j,i)=c(j,i+1)-c(j,i)
1439 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1443 if (itype(i).ne.10) then
1445 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1446 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1454 c----------------------------------------------------------------------------
1455 subroutine read_threadbase
1456 implicit real*8 (a-h,o-z)
1457 include 'DIMENSIONS'
1458 include 'COMMON.IOUNITS'
1459 include 'COMMON.GEO'
1460 include 'COMMON.VAR'
1461 include 'COMMON.INTERACT'
1462 include 'COMMON.LOCAL'
1463 include 'COMMON.NAMES'
1464 include 'COMMON.CHAIN'
1465 include 'COMMON.FFIELD'
1466 include 'COMMON.SBRIDGE'
1467 include 'COMMON.HEADER'
1468 include 'COMMON.CONTROL'
1469 include 'COMMON.DBASE'
1470 include 'COMMON.THREAD'
1471 include 'COMMON.TIME1'
1472 C Read pattern database for threading.
1473 read (icbase,*) nseq
1475 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1476 & nres_base(2,i),nres_base(3,i)
1477 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1479 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1480 c & nres_base(2,i),nres_base(3,i)
1481 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1485 if (weidis.eq.0.0D0) weidis=0.1D0
1494 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1495 write (iout,'(a,i5)') 'nexcl: ',nexcl
1496 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1499 c------------------------------------------------------------------------------
1500 subroutine setup_var
1501 implicit real*8 (a-h,o-z)
1502 include 'DIMENSIONS'
1503 include 'COMMON.IOUNITS'
1504 include 'COMMON.GEO'
1505 include 'COMMON.VAR'
1506 include 'COMMON.INTERACT'
1507 include 'COMMON.LOCAL'
1508 include 'COMMON.NAMES'
1509 include 'COMMON.CHAIN'
1510 include 'COMMON.FFIELD'
1511 include 'COMMON.SBRIDGE'
1512 include 'COMMON.HEADER'
1513 include 'COMMON.CONTROL'
1514 include 'COMMON.DBASE'
1515 include 'COMMON.THREAD'
1516 include 'COMMON.TIME1'
1517 C Set up variable list.
1523 if (itype(i).ne.10) then
1525 ialph(i,1)=nvar+nside
1529 if (indphi.gt.0) then
1531 else if (indback.gt.0) then
1536 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1539 c----------------------------------------------------------------------------
1540 subroutine gen_dist_constr
1541 C Generate CA distance constraints.
1542 implicit real*8 (a-h,o-z)
1543 include 'DIMENSIONS'
1544 include 'COMMON.IOUNITS'
1545 include 'COMMON.GEO'
1546 include 'COMMON.VAR'
1547 include 'COMMON.INTERACT'
1548 include 'COMMON.LOCAL'
1549 include 'COMMON.NAMES'
1550 include 'COMMON.CHAIN'
1551 include 'COMMON.FFIELD'
1552 include 'COMMON.SBRIDGE'
1553 include 'COMMON.HEADER'
1554 include 'COMMON.CONTROL'
1555 include 'COMMON.DBASE'
1556 include 'COMMON.THREAD'
1557 include 'COMMON.TIME1'
1558 dimension itype_pdb(maxres)
1559 common /pizda/ itype_pdb
1561 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1562 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1563 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1565 do i=nstart_sup,nstart_sup+nsup-1
1566 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1567 cd & ' seq_pdb', restyp(itype_pdb(i))
1568 do j=i+2,nstart_sup+nsup-1
1570 ihpb(nhpb)=i+nstart_seq-nstart_sup
1571 jhpb(nhpb)=j+nstart_seq-nstart_sup
1573 dhpb(nhpb)=dist(i,j)
1576 cd write (iout,'(a)') 'Distance constraints:'
1581 cd if (ii.gt.nres) then
1586 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1587 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1588 cd & dhpb(i),forcon(i)
1592 c----------------------------------------------------------------------------
1594 implicit real*8 (a-h,o-z)
1595 include 'DIMENSIONS'
1596 include 'COMMON.MAP'
1597 include 'COMMON.IOUNITS'
1598 character*3 angid(4) /'THE','PHI','ALP','OME'/
1599 character*80 mapcard,ucase
1601 read (inp,'(a)') mapcard
1602 mapcard=ucase(mapcard)
1603 if (index(mapcard,'PHI').gt.0) then
1605 else if (index(mapcard,'THE').gt.0) then
1607 else if (index(mapcard,'ALP').gt.0) then
1609 else if (index(mapcard,'OME').gt.0) then
1612 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1613 stop 'Error - illegal variable spec in MAP card.'
1615 call readi (mapcard,'RES1',res1(imap),0)
1616 call readi (mapcard,'RES2',res2(imap),0)
1617 if (res1(imap).eq.0) then
1618 res1(imap)=res2(imap)
1619 else if (res2(imap).eq.0) then
1620 res2(imap)=res1(imap)
1622 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1624 & 'Error - illegal definition of variable group in MAP.'
1625 stop 'Error - illegal definition of variable group in MAP.'
1627 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1628 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1629 call readi(mapcard,'NSTEP',nstep(imap),0)
1630 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1632 & 'Illegal boundary and/or step size specification in MAP.'
1633 stop 'Illegal boundary and/or step size specification in MAP.'
1638 c----------------------------------------------------------------------------
1639 csa subroutine csaread
1640 csa implicit real*8 (a-h,o-z)
1641 csa include 'DIMENSIONS'
1642 csa include 'COMMON.IOUNITS'
1643 csa include 'COMMON.GEO'
1644 csa include 'COMMON.CSA'
1645 csa include 'COMMON.BANK'
1646 csa include 'COMMON.CONTROL'
1647 csa character*80 ucase
1648 csa character*620 mcmcard
1649 csa call card_concat(mcmcard)
1651 csa call readi(mcmcard,'NCONF',nconf,50)
1652 csa call readi(mcmcard,'NADD',nadd,0)
1653 csa call readi(mcmcard,'JSTART',jstart,1)
1654 csa call readi(mcmcard,'JEND',jend,1)
1655 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1656 csa call readi(mcmcard,'N0',n0,1)
1657 csa call readi(mcmcard,'N1',n1,6)
1658 csa call readi(mcmcard,'N2',n2,4)
1659 csa call readi(mcmcard,'N3',n3,0)
1660 csa call readi(mcmcard,'N4',n4,0)
1661 csa call readi(mcmcard,'N5',n5,0)
1662 csa call readi(mcmcard,'N6',n6,10)
1663 csa call readi(mcmcard,'N7',n7,0)
1664 csa call readi(mcmcard,'N8',n8,0)
1665 csa call readi(mcmcard,'N9',n9,0)
1666 csa call readi(mcmcard,'N14',n14,0)
1667 csa call readi(mcmcard,'N15',n15,0)
1668 csa call readi(mcmcard,'N16',n16,0)
1669 csa call readi(mcmcard,'N17',n17,0)
1670 csa call readi(mcmcard,'N18',n18,0)
1672 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1674 csa call readi(mcmcard,'NDIFF',ndiff,2)
1675 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1676 csa call readi(mcmcard,'IS1',is1,1)
1677 csa call readi(mcmcard,'IS2',is2,8)
1678 csa call readi(mcmcard,'NRAN0',nran0,4)
1679 csa call readi(mcmcard,'NRAN1',nran1,2)
1680 csa call readi(mcmcard,'IRR',irr,1)
1681 csa call readi(mcmcard,'NSEED',nseed,20)
1682 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1683 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1684 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1685 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1686 csa call readi(mcmcard,'ICMAX',icmax,3)
1687 csa call readi(mcmcard,'IRESTART',irestart,0)
1688 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1691 csa call reada(mcmcard,'DELE',dele,20.0d0)
1692 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1693 csa call readi(mcmcard,'IREF',iref,0)
1694 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1695 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1696 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1697 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1698 csa write (iout,*) "NCONF_IN",nconf_in
1701 c----------------------------------------------------------------------------
1702 cfmc subroutine mcmfread
1703 cfmc implicit real*8 (a-h,o-z)
1704 cfmc include 'DIMENSIONS'
1705 cfmc include 'COMMON.MCMF'
1706 cfmc include 'COMMON.IOUNITS'
1707 cfmc include 'COMMON.GEO'
1708 cfmc character*80 ucase
1709 cfmc character*620 mcmcard
1710 cfmc call card_concat(mcmcard)
1712 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1713 cfmc write(iout,*)'MAXRANT=',maxrant
1714 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1715 cfmc write(iout,*)'MAXFAM=',maxfam
1716 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1717 cfmc write(iout,*)'NNET1=',nnet1
1718 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1719 cfmc write(iout,*)'NNET2=',nnet2
1720 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1721 cfmc write(iout,*)'NNET3=',nnet3
1722 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1723 cfmc write(iout,*)'ILASTT=',ilastt
1724 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1725 cfmc write(iout,*)'MAXSTR=',maxstr
1726 cfmc maxstr_f=maxstr/maxfam
1727 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1728 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1729 cfmc write(iout,*)'NMCMF=',nmcmf
1730 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1731 cfmc write(iout,*)'IFOCUS=',ifocus
1732 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1733 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1734 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1735 cfmc write(iout,*)'INTPRT=',intprt
1736 cfmc call readi(mcmcard,'IPRT',iprt,100)
1737 cfmc write(iout,*)'IPRT=',iprt
1738 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1739 cfmc write(iout,*)'IMAXTR=',imaxtr
1740 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1741 cfmc write(iout,*)'MAXEVEN=',maxeven
1742 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1743 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1744 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1745 cfmc write(iout,*)'INIMIN=',inimin
1746 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1747 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1748 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1749 cfmc write(iout,*)'NTHREAD=',nthread
1750 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1751 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1752 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1753 cfmc write(iout,*)'MAXPERT=',maxpert
1754 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1755 cfmc write(iout,*)'IRMSD=',irmsd
1756 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1757 cfmc write(iout,*)'DENEMIN=',denemin
1758 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1759 cfmc write(iout,*)'RCUT1S=',rcut1s
1760 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1761 cfmc write(iout,*)'RCUT1E=',rcut1e
1762 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1763 cfmc write(iout,*)'RCUT2S=',rcut2s
1764 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1765 cfmc write(iout,*)'RCUT2E=',rcut2e
1766 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1767 cfmc write(iout,*)'DPERT1=',d_pert1
1768 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1769 cfmc write(iout,*)'DPERT1A=',d_pert1a
1770 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1771 cfmc write(iout,*)'DPERT2=',d_pert2
1772 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1773 cfmc write(iout,*)'DPERT2A=',d_pert2a
1774 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1775 cfmc write(iout,*)'DPERT2B=',d_pert2b
1776 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1777 cfmc write(iout,*)'DPERT2C=',d_pert2c
1778 cfmc d_pert1=deg2rad*d_pert1
1779 cfmc d_pert1a=deg2rad*d_pert1a
1780 cfmc d_pert2=deg2rad*d_pert2
1781 cfmc d_pert2a=deg2rad*d_pert2a
1782 cfmc d_pert2b=deg2rad*d_pert2b
1783 cfmc d_pert2c=deg2rad*d_pert2c
1784 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1785 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1786 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1787 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1788 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1789 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1790 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1791 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1792 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1793 cfmc write(iout,*)'RCUTINI=',rcutini
1794 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1795 cfmc write(iout,*)'GRAT=',grat
1796 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1797 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1801 c----------------------------------------------------------------------------
1803 implicit real*8 (a-h,o-z)
1804 include 'DIMENSIONS'
1805 include 'COMMON.MCM'
1806 include 'COMMON.MCE'
1807 include 'COMMON.IOUNITS'
1809 character*320 mcmcard
1810 call card_concat(mcmcard)
1811 call readi(mcmcard,'MAXACC',maxacc,100)
1812 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1813 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1814 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1815 call readi(mcmcard,'MAXREPM',maxrepm,200)
1816 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1817 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1818 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1819 call reada(mcmcard,'E_UP',e_up,5.0D0)
1820 call reada(mcmcard,'DELTE',delte,0.1D0)
1821 call readi(mcmcard,'NSWEEP',nsweep,5)
1822 call readi(mcmcard,'NSTEPH',nsteph,0)
1823 call readi(mcmcard,'NSTEPC',nstepc,0)
1824 call reada(mcmcard,'TMIN',tmin,298.0D0)
1825 call reada(mcmcard,'TMAX',tmax,298.0D0)
1826 call readi(mcmcard,'NWINDOW',nwindow,0)
1827 call readi(mcmcard,'PRINT_MC',print_mc,0)
1828 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1829 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1830 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1831 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1832 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1833 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1834 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1835 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1836 if (nwindow.gt.0) then
1837 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1839 winlen(i)=winend(i)-winstart(i)+1
1842 if (tmax.lt.tmin) tmax=tmin
1843 if (tmax.eq.tmin) then
1847 if (nstepc.gt.0 .and. nsteph.gt.0) then
1848 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1849 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1851 C Probabilities of different move types
1852 sumpro_type(0)=0.0D0
1853 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1854 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1855 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1856 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1857 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1858 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1859 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1861 print *,'i',i,' sumprotype',sumpro_type(i)
1862 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1863 print *,'i',i,' sumprotype',sumpro_type(i)
1867 c----------------------------------------------------------------------------
1868 subroutine read_minim
1869 implicit real*8 (a-h,o-z)
1870 include 'DIMENSIONS'
1871 include 'COMMON.MINIM'
1872 include 'COMMON.IOUNITS'
1874 character*320 minimcard
1875 call card_concat(minimcard)
1876 call readi(minimcard,'MAXMIN',maxmin,2000)
1877 call readi(minimcard,'MAXFUN',maxfun,5000)
1878 call readi(minimcard,'MINMIN',minmin,maxmin)
1879 call readi(minimcard,'MINFUN',minfun,maxmin)
1880 call reada(minimcard,'TOLF',tolf,1.0D-2)
1881 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1882 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1883 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1884 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1885 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1886 & 'Options in energy minimization:'
1887 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1888 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1889 & 'MinMin:',MinMin,' MinFun:',MinFun,
1890 & ' TolF:',TolF,' RTolF:',RTolF
1893 c----------------------------------------------------------------------------
1894 subroutine read_angles(kanal,*)
1895 implicit real*8 (a-h,o-z)
1896 include 'DIMENSIONS'
1897 include 'COMMON.GEO'
1898 include 'COMMON.VAR'
1899 include 'COMMON.CHAIN'
1900 include 'COMMON.IOUNITS'
1901 include 'COMMON.CONTROL'
1902 c Read angles from input
1904 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1905 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1906 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1907 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1910 c 9/7/01 avoid 180 deg valence angle
1911 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1913 theta(i)=deg2rad*theta(i)
1914 phi(i)=deg2rad*phi(i)
1915 alph(i)=deg2rad*alph(i)
1916 omeg(i)=deg2rad*omeg(i)
1921 c----------------------------------------------------------------------------
1922 subroutine reada(rekord,lancuch,wartosc,default)
1924 character*(*) rekord,lancuch
1925 double precision wartosc,default
1928 iread=index(rekord,lancuch)
1929 if (iread.eq.0) then
1933 iread=iread+ilen(lancuch)+1
1934 read (rekord(iread:),*,err=10,end=10) wartosc
1939 c----------------------------------------------------------------------------
1940 subroutine readi(rekord,lancuch,wartosc,default)
1942 character*(*) rekord,lancuch
1943 integer wartosc,default
1946 iread=index(rekord,lancuch)
1947 if (iread.eq.0) then
1951 iread=iread+ilen(lancuch)+1
1952 read (rekord(iread:),*,err=10,end=10) wartosc
1957 c----------------------------------------------------------------------------
1958 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1961 integer tablica(dim),default
1962 character*(*) rekord,lancuch
1969 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1970 if (iread.eq.0) return
1971 iread=iread+ilen(lancuch)+1
1972 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1975 c----------------------------------------------------------------------------
1976 subroutine multreada(rekord,lancuch,tablica,dim,default)
1979 double precision tablica(dim),default
1980 character*(*) rekord,lancuch
1987 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1988 if (iread.eq.0) return
1989 iread=iread+ilen(lancuch)+1
1990 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1993 c----------------------------------------------------------------------------
1994 subroutine openunits
1995 implicit real*8 (a-h,o-z)
1996 include 'DIMENSIONS'
1999 character*16 form,nodename
2002 include 'COMMON.SETUP'
2003 include 'COMMON.IOUNITS'
2005 include 'COMMON.CONTROL'
2006 integer lenpre,lenpot,ilen,lentmp
2008 character*3 out1file_text,ucase
2011 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2012 call getenv_loc("PREFIX",prefix)
2014 call getenv_loc("POT",pot)
2015 call getenv_loc("DIRTMP",tmpdir)
2016 call getenv_loc("CURDIR",curdir)
2017 call getenv_loc("OUT1FILE",out1file_text)
2018 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2019 out1file_text=ucase(out1file_text)
2020 if (out1file_text(1:1).eq."Y") then
2023 out1file=fg_rank.gt.0
2028 if (lentmp.gt.0) then
2029 write (*,'(80(1h!))')
2030 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2031 write (*,'(80(1h!))')
2032 write (*,*)"All output files will be on node /tmp directory."
2034 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2035 if (me.eq.king) then
2036 write (*,*) "The master node is ",nodename
2037 else if (fg_rank.eq.0) then
2038 write (*,*) "I am the CG slave node ",nodename
2040 write (*,*) "I am the FG slave node ",nodename
2043 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2044 lenpre = lentmp+lenpre+1
2046 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2047 C Get the names and open the input files
2048 #if defined(WINIFL) || defined(WINPGI)
2049 open(1,file=pref_orig(:ilen(pref_orig))//
2050 & '.inp',status='old',readonly,shared)
2051 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2052 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2053 C Get parameter filenames and open the parameter files.
2054 call getenv_loc('BONDPAR',bondname)
2055 open (ibond,file=bondname,status='old',readonly,shared)
2056 call getenv_loc('THETPAR',thetname)
2057 open (ithep,file=thetname,status='old',readonly,shared)
2059 call getenv_loc('THETPARPDB',thetname_pdb)
2060 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2062 call getenv_loc('ROTPAR',rotname)
2063 open (irotam,file=rotname,status='old',readonly,shared)
2065 call getenv_loc('ROTPARPDB',rotname_pdb)
2066 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2068 call getenv_loc('TORPAR',torname)
2069 open (itorp,file=torname,status='old',readonly,shared)
2070 call getenv_loc('TORDPAR',tordname)
2071 open (itordp,file=tordname,status='old',readonly,shared)
2072 call getenv_loc('FOURIER',fouriername)
2073 open (ifourier,file=fouriername,status='old',readonly,shared)
2074 call getenv_loc('ELEPAR',elename)
2075 open (ielep,file=elename,status='old',readonly,shared)
2076 call getenv_loc('SIDEPAR',sidename)
2077 open (isidep,file=sidename,status='old',readonly,shared)
2078 #elif (defined CRAY) || (defined AIX)
2079 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2081 c print *,"Processor",myrank," opened file 1"
2082 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2083 c print *,"Processor",myrank," opened file 9"
2084 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2085 C Get parameter filenames and open the parameter files.
2086 call getenv_loc('BONDPAR',bondname)
2087 open (ibond,file=bondname,status='old',action='read')
2088 c print *,"Processor",myrank," opened file IBOND"
2089 call getenv_loc('THETPAR',thetname)
2090 open (ithep,file=thetname,status='old',action='read')
2091 c print *,"Processor",myrank," opened file ITHEP"
2093 call getenv_loc('THETPARPDB',thetname_pdb)
2094 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2096 call getenv_loc('ROTPAR',rotname)
2097 open (irotam,file=rotname,status='old',action='read')
2098 c print *,"Processor",myrank," opened file IROTAM"
2100 call getenv_loc('ROTPARPDB',rotname_pdb)
2101 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2103 call getenv_loc('TORPAR',torname)
2104 open (itorp,file=torname,status='old',action='read')
2105 c print *,"Processor",myrank," opened file ITORP"
2106 call getenv_loc('TORDPAR',tordname)
2107 open (itordp,file=tordname,status='old',action='read')
2108 c print *,"Processor",myrank," opened file ITORDP"
2109 call getenv_loc('SCCORPAR',sccorname)
2110 open (isccor,file=sccorname,status='old',action='read')
2111 c print *,"Processor",myrank," opened file ISCCOR"
2112 call getenv_loc('FOURIER',fouriername)
2113 open (ifourier,file=fouriername,status='old',action='read')
2114 c print *,"Processor",myrank," opened file IFOURIER"
2115 call getenv_loc('ELEPAR',elename)
2116 open (ielep,file=elename,status='old',action='read')
2117 c print *,"Processor",myrank," opened file IELEP"
2118 call getenv_loc('SIDEPAR',sidename)
2119 open (isidep,file=sidename,status='old',action='read')
2120 c print *,"Processor",myrank," opened file ISIDEP"
2121 c print *,"Processor",myrank," opened parameter files"
2123 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2124 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2125 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2126 C Get parameter filenames and open the parameter files.
2127 call getenv_loc('BONDPAR',bondname)
2128 open (ibond,file=bondname,status='old')
2129 call getenv_loc('THETPAR',thetname)
2130 open (ithep,file=thetname,status='old')
2132 call getenv_loc('THETPARPDB',thetname_pdb)
2133 open (ithep_pdb,file=thetname_pdb,status='old')
2135 call getenv_loc('ROTPAR',rotname)
2136 open (irotam,file=rotname,status='old')
2138 call getenv_loc('ROTPARPDB',rotname_pdb)
2139 open (irotam_pdb,file=rotname_pdb,status='old')
2141 call getenv_loc('TORPAR',torname)
2142 open (itorp,file=torname,status='old')
2143 call getenv_loc('TORDPAR',tordname)
2144 open (itordp,file=tordname,status='old')
2145 call getenv_loc('SCCORPAR',sccorname)
2146 open (isccor,file=sccorname,status='old')
2147 call getenv_loc('FOURIER',fouriername)
2148 open (ifourier,file=fouriername,status='old')
2149 call getenv_loc('ELEPAR',elename)
2150 open (ielep,file=elename,status='old')
2151 call getenv_loc('SIDEPAR',sidename)
2152 open (isidep,file=sidename,status='old')
2154 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2156 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2157 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2158 C Get parameter filenames and open the parameter files.
2159 call getenv_loc('BONDPAR',bondname)
2160 open (ibond,file=bondname,status='old',action='read')
2161 call getenv_loc('THETPAR',thetname)
2162 open (ithep,file=thetname,status='old',action='read')
2164 call getenv_loc('THETPARPDB',thetname_pdb)
2165 print *,"thetname_pdb ",thetname_pdb
2166 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2167 print *,ithep_pdb," opened"
2169 call getenv_loc('ROTPAR',rotname)
2170 open (irotam,file=rotname,status='old',action='read')
2172 call getenv_loc('ROTPARPDB',rotname_pdb)
2173 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2175 call getenv_loc('TORPAR',torname)
2176 open (itorp,file=torname,status='old',action='read')
2177 call getenv_loc('TORDPAR',tordname)
2178 open (itordp,file=tordname,status='old',action='read')
2179 call getenv_loc('SCCORPAR',sccorname)
2180 open (isccor,file=sccorname,status='old',action='read')
2181 call getenv_loc('FOURIER',fouriername)
2182 open (ifourier,file=fouriername,status='old',action='read')
2183 call getenv_loc('ELEPAR',elename)
2184 open (ielep,file=elename,status='old',action='read')
2185 call getenv_loc('SIDEPAR',sidename)
2186 open (isidep,file=sidename,status='old',action='read')
2190 C 8/9/01 In the newest version SCp interaction constants are read from a file
2191 C Use -DOLDSCP to use hard-coded constants instead.
2193 call getenv_loc('SCPPAR',scpname)
2194 #if defined(WINIFL) || defined(WINPGI)
2195 open (iscpp,file=scpname,status='old',readonly,shared)
2196 #elif (defined CRAY) || (defined AIX)
2197 open (iscpp,file=scpname,status='old',action='read')
2199 open (iscpp,file=scpname,status='old')
2201 open (iscpp,file=scpname,status='old',action='read')
2204 call getenv_loc('PATTERN',patname)
2205 #if defined(WINIFL) || defined(WINPGI)
2206 open (icbase,file=patname,status='old',readonly,shared)
2207 #elif (defined CRAY) || (defined AIX)
2208 open (icbase,file=patname,status='old',action='read')
2210 open (icbase,file=patname,status='old')
2212 open (icbase,file=patname,status='old',action='read')
2215 C Open output file only for CG processes
2216 c print *,"Processor",myrank," fg_rank",fg_rank
2217 if (fg_rank.eq.0) then
2219 if (nodes.eq.1) then
2222 npos = dlog10(dfloat(nodes-1))+1
2224 if (npos.lt.3) npos=3
2225 write (liczba,'(i1)') npos
2226 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2228 write (liczba,form) me
2229 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2230 & liczba(:ilen(liczba))
2231 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2233 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2235 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2236 & liczba(:ilen(liczba))//'.mol2'
2237 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2238 & liczba(:ilen(liczba))//'.stat'
2240 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2241 & //liczba(:ilen(liczba))//'.stat')
2242 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2245 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2246 & liczba(:ilen(liczba))//'.const'
2251 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2252 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2253 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2254 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2255 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2257 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2259 rest2name=prefix(:ilen(prefix))//'.rst'
2261 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2264 #if defined(AIX) || defined(PGI)
2265 if (me.eq.king .or. .not. out1file)
2266 & open(iout,file=outname,status='unknown')
2269 if (fg_rank.gt.0) then
2270 write (liczba,'(i3.3)') myrank/nfgtasks
2271 write (ll,'(bz,i3.3)') fg_rank
2272 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2278 open(igeom,file=intname,status='unknown',position='append')
2279 open(ipdb,file=pdbname,status='unknown')
2280 open(imol2,file=mol2name,status='unknown')
2281 open(istat,file=statname,status='unknown',position='append')
2283 c1out open(iout,file=outname,status='unknown')
2286 if (me.eq.king .or. .not.out1file)
2287 & open(iout,file=outname,status='unknown')
2290 if (fg_rank.gt.0) then
2291 print "Processor",fg_rank," opening output file"
2292 write (liczba,'(i3.3)') myrank/nfgtasks
2293 write (ll,'(bz,i3.3)') fg_rank
2294 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2300 open(igeom,file=intname,status='unknown',access='append')
2301 open(ipdb,file=pdbname,status='unknown')
2302 open(imol2,file=mol2name,status='unknown')
2303 open(istat,file=statname,status='unknown',access='append')
2305 c1out open(iout,file=outname,status='unknown')
2308 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2309 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2310 csa csa_history=prefix(:lenpre)//'.CSA.history'
2311 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2312 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2313 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2314 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2315 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2316 csa csa_int=prefix(:lenpre)//'.int'
2317 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2318 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2319 csa csa_in=prefix(:lenpre)//'.CSA.in'
2320 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2323 write (iout,'(80(1h-))')
2324 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2325 write (iout,'(80(1h-))')
2326 write (iout,*) "Input file : ",
2327 & pref_orig(:ilen(pref_orig))//'.inp'
2328 write (iout,*) "Output file : ",
2329 & outname(:ilen(outname))
2331 write (iout,*) "Sidechain potential file : ",
2332 & sidename(:ilen(sidename))
2334 write (iout,*) "SCp potential file : ",
2335 & scpname(:ilen(scpname))
2337 write (iout,*) "Electrostatic potential file : ",
2338 & elename(:ilen(elename))
2339 write (iout,*) "Cumulant coefficient file : ",
2340 & fouriername(:ilen(fouriername))
2341 write (iout,*) "Torsional parameter file : ",
2342 & torname(:ilen(torname))
2343 write (iout,*) "Double torsional parameter file : ",
2344 & tordname(:ilen(tordname))
2345 write (iout,*) "SCCOR parameter file : ",
2346 & sccorname(:ilen(sccorname))
2347 write (iout,*) "Bond & inertia constant file : ",
2348 & bondname(:ilen(bondname))
2349 write (iout,*) "Bending parameter file : ",
2350 & thetname(:ilen(thetname))
2351 write (iout,*) "Rotamer parameter file : ",
2352 & rotname(:ilen(rotname))
2353 write (iout,*) "Threading database : ",
2354 & patname(:ilen(patname))
2356 &write (iout,*)" DIRTMP : ",
2358 write (iout,'(80(1h-))')
2362 c----------------------------------------------------------------------------
2363 subroutine card_concat(card)
2364 implicit real*8 (a-h,o-z)
2365 include 'DIMENSIONS'
2366 include 'COMMON.IOUNITS'
2368 character*80 karta,ucase
2370 read (inp,'(a)') karta
2373 do while (karta(80:80).eq.'&')
2374 card=card(:ilen(card)+1)//karta(:79)
2375 read (inp,'(a)') karta
2378 card=card(:ilen(card)+1)//karta
2381 c----------------------------------------------------------------------------------
2383 implicit real*8 (a-h,o-z)
2384 include 'DIMENSIONS'
2385 include 'COMMON.CHAIN'
2386 include 'COMMON.IOUNITS'
2388 open(irest2,file=rest2name,status='unknown')
2389 read(irest2,*) totT,EK,potE,totE,t_bath
2391 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2394 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2397 read (irest2,*) iset
2402 c---------------------------------------------------------------------------------
2403 subroutine read_fragments
2404 implicit real*8 (a-h,o-z)
2405 include 'DIMENSIONS'
2409 include 'COMMON.SETUP'
2410 include 'COMMON.CHAIN'
2411 include 'COMMON.IOUNITS'
2413 include 'COMMON.CONTROL'
2414 read(inp,*) nset,nfrag,npair,nfrag_back
2415 if(me.eq.king.or..not.out1file)
2416 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2417 & " nfrag_back",nfrag_back
2419 read(inp,*) mset(iset)
2421 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2423 if(me.eq.king.or..not.out1file)
2424 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2425 & ifrag(2,i,iset), qinfrag(i,iset)
2428 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2430 if(me.eq.king.or..not.out1file)
2431 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2432 & ipair(2,i,iset), qinpair(i,iset)
2435 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2436 & wfrag_back(3,i,iset),
2437 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2438 if(me.eq.king.or..not.out1file)
2439 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2440 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2445 c-------------------------------------------------------------------------------
2446 subroutine read_dist_constr
2447 implicit real*8 (a-h,o-z)
2448 include 'DIMENSIONS'
2452 include 'COMMON.SETUP'
2453 include 'COMMON.CONTROL'
2454 include 'COMMON.CHAIN'
2455 include 'COMMON.IOUNITS'
2456 include 'COMMON.SBRIDGE'
2457 integer ifrag_(2,100),ipair_(2,100)
2458 double precision wfrag_(100),wpair_(100)
2459 character*500 controlcard
2460 c write (iout,*) "Calling read_dist_constr"
2461 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2463 call card_concat(controlcard)
2464 call readi(controlcard,"NFRAG",nfrag_,0)
2465 call readi(controlcard,"NPAIR",npair_,0)
2466 call readi(controlcard,"NDIST",ndist_,0)
2467 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2468 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2469 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2470 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2471 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2472 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2473 c write (iout,*) "IFRAG"
2475 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2477 c write (iout,*) "IPAIR"
2479 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2481 if (.not.refstr .and. nfrag.gt.0) then
2483 & "ERROR: no reference structure to compute distance restraints"
2485 & "Restraints must be specified explicitly (NDIST=number)"
2488 if (nfrag.lt.2 .and. npair.gt.0) then
2489 write (iout,*) "ERROR: Less than 2 fragments specified",
2490 & " but distance restraints between pairs requested"
2495 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2496 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2497 & ifrag_(2,i)=nstart_sup+nsup-1
2498 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2500 if (wfrag_(i).gt.0.0d0) then
2501 do j=ifrag_(1,i),ifrag_(2,i)-1
2502 do k=j+1,ifrag_(2,i)
2503 c write (iout,*) "j",j," k",k
2505 if (constr_dist.eq.1) then
2510 forcon(nhpb)=wfrag_(i)
2511 else if (constr_dist.eq.2) then
2512 if (ddjk.le.dist_cut) then
2517 forcon(nhpb)=wfrag_(i)
2524 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2527 if (.not.out1file .or. me.eq.king)
2528 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2529 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2531 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2532 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2539 if (wpair_(i).gt.0.0d0) then
2547 do j=ifrag_(1,ii),ifrag_(2,ii)
2548 do k=ifrag_(1,jj),ifrag_(2,jj)
2552 forcon(nhpb)=wpair_(i)
2553 dhpb(nhpb)=dist(j,k)
2555 if (.not.out1file .or. me.eq.king)
2556 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2557 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2559 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2560 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2567 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2568 & ibecarb(i),forcon(nhpb+1)
2569 if (forcon(nhpb+1).gt.0.0d0) then
2571 if (ibecarb(i).gt.0) then
2572 ihpb(i)=ihpb(i)+nres
2573 jhpb(i)=jhpb(i)+nres
2575 if (dhpb(nhpb).eq.0.0d0)
2576 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2580 if (.not.out1file .or. me.eq.king) then
2583 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2584 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2592 c-------------------------------------------------------------------------------
2594 subroutine flush(iu)
2599 subroutine flush(iu)
2604 c------------------------------------------------------------------------------
2605 subroutine copy_to_tmp(source)
2606 include "DIMENSIONS"
2607 include "COMMON.IOUNITS"
2608 character*(*) source
2609 character* 256 tmpfile
2613 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2614 inquire(file=tmpfile,exist=ex)
2616 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2617 & " to temporary directory..."
2618 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2619 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2623 c------------------------------------------------------------------------------
2624 subroutine move_from_tmp(source)
2625 include "DIMENSIONS"
2626 include "COMMON.IOUNITS"
2627 character*(*) source
2630 write (*,*) "Moving ",source(:ilen(source)),
2631 & " from temporary directory to working directory"
2632 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2633 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2636 c------------------------------------------------------------------------------
2637 subroutine random_init(seed)
2639 C Initialize random number generator
2641 implicit real*8 (a-h,o-z)
2642 include 'DIMENSIONS'
2648 logical OKRandom, prng_restart
2650 integer iseed_array(4)
2652 include 'COMMON.IOUNITS'
2653 include 'COMMON.TIME1'
2654 include 'COMMON.THREAD'
2655 include 'COMMON.SBRIDGE'
2656 include 'COMMON.CONTROL'
2657 include 'COMMON.MCM'
2658 include 'COMMON.MAP'
2659 include 'COMMON.HEADER'
2660 csa include 'COMMON.CSA'
2661 include 'COMMON.CHAIN'
2662 include 'COMMON.MUCA'
2664 include 'COMMON.FFIELD'
2665 include 'COMMON.SETUP'
2666 iseed=-dint(dabs(seed))
2667 if (iseed.eq.0) then
2668 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2669 & 'Random seed undefined. The program will stop.'
2670 write (*,'(/80(1h*)/20x,a/80(1h*))')
2671 & 'Random seed undefined. The program will stop.'
2673 call mpi_finalize(mpi_comm_world,ierr)
2675 stop 'Bad random seed.'
2678 if (fg_rank.eq.0) then
2682 if(me.eq.king .or. .not. out1file)
2683 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2684 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2685 OKRandom = prng_restart(me,iseedi8)
2688 tmp=65536.0d0**(4-i)
2689 iseed_array(i) = dint(seed/tmp)
2690 seed=seed-iseed_array(i)*tmp
2692 if(me.eq.king .or. .not. out1file)
2693 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2694 & (iseed_array(i),i=1,4)
2695 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2696 & (iseed_array(i),i=1,4)
2697 OKRandom = prng_restart(me,iseed_array)
2700 r1=ran_number(0.0D0,1.0D0)
2701 if(me.eq.king .or. .not. out1file)
2702 & write (iout,*) 'ran_num',r1
2703 if (r1.lt.0.0d0) OKRandom=.false.
2705 if (.not.OKRandom) then
2706 write (iout,*) 'PRNG IS NOT WORKING!!!'
2707 print *,'PRNG IS NOT WORKING!!!'
2710 call mpi_abort(mpi_comm_world,error_msg,ierr)
2713 write (iout,*) 'too many processors for parallel prng'
2714 write (*,*) 'too many processors for parallel prng'
2722 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)