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 call reada(weightcard,"DTRISS",dtriss,1.0D0)
864 call reada(weightcard,"ATRISS",atriss,0.3D0)
865 call reada(weightcard,"BTRISS",btriss,0.02D0)
866 call reada(weightcard,"CTRISS",ctriss,1.0D0)
867 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
869 dyn_ss_mask(i)=.false.
873 dyn_ssbond_ij(i,j)=1.0d300
876 call reada(weightcard,"HT",Ht,0.0D0)
878 ss_depth=ebr/wsc-0.25*eps(1,1)
879 Ht=Ht/wsc-0.25*eps(1,1)
880 akcm=akcm*wstrain/wsc
881 akth=akth*wstrain/wsc
882 akct=akct*wstrain/wsc
883 v1ss=v1ss*wstrain/wsc
884 v2ss=v2ss*wstrain/wsc
885 v3ss=v3ss*wstrain/wsc
887 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
890 if(me.eq.king.or..not.out1file) then
891 write (iout,*) "Parameters of the SS-bond potential:"
892 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
894 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
895 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
896 write (iout,*)" HT",Ht
897 print *,'indpdb=',indpdb,' pdbref=',pdbref
899 if (indpdb.gt.0 .or. pdbref) then
900 read(inp,'(a)') pdbfile
901 if(me.eq.king.or..not.out1file)
902 & write (iout,'(2a)') 'PDB data will be read from file ',
903 & pdbfile(:ilen(pdbfile))
904 open(ipdbin,file=pdbfile,status='old',err=33)
906 33 write (iout,'(a)') 'Error opening PDB file.'
909 c print *,'Begin reading pdb data'
911 c print *,'Finished reading pdb data'
912 if(me.eq.king.or..not.out1file)
913 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
914 & ' nstart_sup=',nstart_sup
916 itype_pdb(i)=itype(i)
920 nct=nstart_sup+nsup-1
921 call contact(.false.,ncont_ref,icont_ref,co)
924 C Following 2 lines for diagnostics; comment out if not needed
925 write (iout,*) "Before sideadd"
927 if(me.eq.king.or..not.out1file)
928 & write(iout,*)'Adding sidechains'
935 do while (fail.and.nsi.le.maxsi)
936 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
939 if(fail) write(iout,*)'Adding sidechain failed for res ',
940 & i,' after ',nsi,' trials'
943 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
946 C Following 2 lines for diagnostics; comment out if not needed
947 c write (iout,*) "After sideadd"
950 if (indpdb.eq.0) then
951 C Read sequence if not taken from the pdb file.
953 c print *,'nres=',nres
954 if (iscode.gt.0) then
955 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
957 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
959 C Convert sequence to numeric code
961 itype(i)=rescode(i,sequence(i),iscode)
963 C Assign initial virtual bond lengths
969 vbld(i+nres)=dsc(itype(i))
970 vbld_inv(i+nres)=dsc_inv(itype(i))
971 c write (iout,*) "i",i," itype",itype(i),
972 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
976 c print '(20i4)',(itype(i),i=1,nres)
979 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
981 if (itype(i).eq.ntyp1) then
985 else if (itype(i+1).ne.20) then
987 else if (itype(i).ne.20) then
994 if(me.eq.king.or..not.out1file)then
995 write (iout,*) "ITEL"
997 write (iout,*) i,itype(i),itel(i)
999 print *,'Call Read_Bridge.'
1002 C 8/13/98 Set limits to generating the dihedral angles
1007 read (inp,*) ndih_constr
1008 if (ndih_constr.gt.0) then
1010 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1011 if(me.eq.king.or..not.out1file)then
1013 & 'There are',ndih_constr,' constraints on phi angles.'
1015 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1019 phi0(i)=deg2rad*phi0(i)
1020 drange(i)=deg2rad*drange(i)
1022 if(me.eq.king.or..not.out1file)
1023 & write (iout,*) 'FTORS',ftors
1026 phibound(1,ii) = phi0(i)-drange(i)
1027 phibound(2,ii) = phi0(i)+drange(i)
1032 if (me.eq.king) then
1034 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1036 write (iout,'(a3,i5,2f10.1)')
1037 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1043 cd print *,'NNT=',NNT,' NCT=',NCT
1044 if (itype(1).eq.ntyp1) nnt=2
1045 if (itype(nres).eq.ntyp1) nct=nct-1
1047 if(me.eq.king.or..not.out1file)
1048 & write (iout,'(a,i3)') 'nsup=',nsup
1050 if (nsup.le.(nct-nnt+1)) then
1051 do i=0,nct-nnt+1-nsup
1052 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1058 & 'Error - sequences to be superposed do not match.'
1061 do i=0,nsup-(nct-nnt+1)
1062 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1064 nstart_sup=nstart_sup+i
1070 & 'Error - sequences to be superposed do not match.'
1073 if (nsup.eq.0) nsup=nct-nnt
1074 if (nstart_sup.eq.0) nstart_sup=nnt
1075 if (nstart_seq.eq.0) nstart_seq=nnt
1076 if(me.eq.king.or..not.out1file)
1077 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1078 & ' nstart_seq=',nstart_seq
1080 c--- Zscore rms -------
1081 if (nz_start.eq.0) nz_start=nnt
1082 if (nz_end.eq.0 .and. nsup.gt.0) then
1084 else if (nz_end.eq.0) then
1087 if(me.eq.king.or..not.out1file)then
1088 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1089 write (iout,*) 'IZ_SC=',iz_sc
1091 c----------------------
1094 if (.not.pdbref) then
1095 call read_angles(inp,*38)
1097 38 write (iout,'(a)') 'Error reading reference structure.'
1099 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1100 stop 'Error reading reference structure'
1104 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1113 call contact(.true.,ncont_ref,icont_ref,co)
1115 if(me.eq.king.or..not.out1file)
1116 & write (iout,*) 'Contact order:',co
1118 if(me.eq.king.or..not.out1file)
1119 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1122 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1124 if(me.eq.king.or..not.out1file)
1125 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1126 & icont_ref(1,i),' ',
1127 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1131 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1132 if (constr_dist.gt.0) then
1133 call read_dist_constr
1135 if (nhpb.gt.0) call hpb_partition
1136 c write (iout,*) "After read_dist_constr nhpb",nhpb
1138 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1139 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1140 & modecalc.ne.10) then
1141 C If input structure hasn't been supplied from the PDB file read or generate
1143 if (iranconf.eq.0 .and. .not. extconf) then
1144 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1145 & write (iout,'(a)') 'Initial geometry will be read in.'
1147 read(inp,'(8f10.5)',end=36,err=36)
1148 & ((c(l,k),l=1,3),k=1,nres),
1149 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1150 call int_from_cart1(.false.)
1153 dc(j,i)=c(j,i+1)-c(j,i)
1154 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1158 if (itype(i).ne.10) then
1160 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1161 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1167 call read_angles(inp,*36)
1170 36 write (iout,'(a)') 'Error reading angle file.'
1172 call mpi_finalize( MPI_COMM_WORLD,IERR )
1174 stop 'Error reading angle file.'
1176 else if (extconf) then
1177 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1178 & write (iout,'(a)') 'Extended chain initial geometry.'
1180 theta(i)=90d0*deg2rad
1183 phi(i)=180d0*deg2rad
1186 alph(i)=110d0*deg2rad
1189 omeg(i)=-120d0*deg2rad
1193 if(me.eq.king.or..not.out1file)
1194 & write (iout,'(a)') 'Random-generated initial geometry.'
1198 if (me.eq.king .or. fg_rank.eq.0 .and. (
1199 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1203 call gen_rand_conf(itmp,*30)
1205 30 write (iout,*) 'Failed to generate random conformation',
1206 & ', itrial=',itrial
1207 write (*,*) 'Processor:',me,
1208 & ' Failed to generate random conformation',
1218 write (iout,'(a,i3,a)') 'Processor:',me,
1219 & ' error in generating random conformation.'
1220 write (*,'(a,i3,a)') 'Processor:',me,
1221 & ' error in generating random conformation.'
1224 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1231 elseif (modecalc.eq.4) then
1232 read (inp,'(a)') intinname
1233 open (intin,file=intinname,status='old',err=333)
1234 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1235 & write (iout,'(a)') 'intinname',intinname
1236 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1238 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1240 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1242 stop 'Error opening angle file.'
1246 C Generate distance constraints, if the PDB structure is to be regularized.
1247 if (nthread.gt.0) then
1248 call read_threadbase
1251 if (me.eq.king .or. .not. out1file)
1253 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1254 write (iout,'(/a,i3,a)')
1255 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1256 write (iout,'(20i4)') (iss(i),i=1,ns)
1258 write(iout,*)"Running with dynamic disulfide-bond formation"
1260 write (iout,'(/a/)') 'Pre-formed links are:'
1266 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1267 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1273 if (ns.gt.0.and.dyn_ss) then
1277 forcon(i-nss)=forcon(i)
1284 dyn_ss_mask(iss(i))=.true.
1287 if (i2ndstr.gt.0) call secstrp2dihc
1288 c call geom_to_var(nvar,x)
1289 c call etotal(energia(0))
1290 c call enerprint(energia(0))
1291 c call briefout(0,etot)
1293 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1294 cd write (iout,'(a)') 'Variable list:'
1295 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1297 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1298 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1299 & 'Processor',myrank,': end reading molecular data.'
1303 c--------------------------------------------------------------------------
1304 logical function seq_comp(itypea,itypeb,length)
1306 integer length,itypea(length),itypeb(length)
1309 if (itypea(i).ne.itypeb(i)) then
1317 c-----------------------------------------------------------------------------
1318 subroutine read_bridge
1319 C Read information about disulfide bridges.
1320 implicit real*8 (a-h,o-z)
1321 include 'DIMENSIONS'
1325 include 'COMMON.IOUNITS'
1326 include 'COMMON.GEO'
1327 include 'COMMON.VAR'
1328 include 'COMMON.INTERACT'
1329 include 'COMMON.LOCAL'
1330 include 'COMMON.NAMES'
1331 include 'COMMON.CHAIN'
1332 include 'COMMON.FFIELD'
1333 include 'COMMON.SBRIDGE'
1334 include 'COMMON.HEADER'
1335 include 'COMMON.CONTROL'
1336 include 'COMMON.DBASE'
1337 include 'COMMON.THREAD'
1338 include 'COMMON.TIME1'
1339 include 'COMMON.SETUP'
1340 C Read bridging residues.
1341 read (inp,*) ns,(iss(i),i=1,ns)
1343 if(me.eq.king.or..not.out1file)
1344 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1345 C Check whether the specified bridging residues are cystines.
1347 if (itype(iss(i)).ne.1) then
1348 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1349 & 'Do you REALLY think that the residue ',
1350 & restyp(itype(iss(i))),i,
1351 & ' can form a disulfide bridge?!!!'
1352 write (*,'(2a,i3,a)')
1353 & 'Do you REALLY think that the residue ',
1354 & restyp(itype(iss(i))),i,
1355 & ' can form a disulfide bridge?!!!'
1357 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1362 C Read preformed bridges.
1364 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1366 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1369 C Check if the residues involved in bridges are in the specified list of
1370 C bridging residues.
1373 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1374 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1375 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1376 & ' contains residues present in other pairs.'
1377 write (*,'(a,i3,a)') 'Disulfide pair',i,
1378 & ' contains residues present in other pairs.'
1380 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1386 if (ihpb(i).eq.iss(j)) goto 10
1388 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1391 if (jhpb(i).eq.iss(j)) goto 20
1393 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1399 ihpb(i)=ihpb(i)+nres
1400 jhpb(i)=jhpb(i)+nres
1406 c----------------------------------------------------------------------------
1407 subroutine read_x(kanal,*)
1408 implicit real*8 (a-h,o-z)
1409 include 'DIMENSIONS'
1410 include 'COMMON.GEO'
1411 include 'COMMON.VAR'
1412 include 'COMMON.CHAIN'
1413 include 'COMMON.IOUNITS'
1414 include 'COMMON.CONTROL'
1415 include 'COMMON.LOCAL'
1416 include 'COMMON.INTERACT'
1417 c Read coordinates from input
1419 read(kanal,'(8f10.5)',end=10,err=10)
1420 & ((c(l,k),l=1,3),k=1,nres),
1421 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1424 c(j,2*nres)=c(j,nres)
1426 call int_from_cart1(.false.)
1429 dc(j,i)=c(j,i+1)-c(j,i)
1430 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1434 if (itype(i).ne.10) then
1436 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1437 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1445 c----------------------------------------------------------------------------
1446 subroutine read_threadbase
1447 implicit real*8 (a-h,o-z)
1448 include 'DIMENSIONS'
1449 include 'COMMON.IOUNITS'
1450 include 'COMMON.GEO'
1451 include 'COMMON.VAR'
1452 include 'COMMON.INTERACT'
1453 include 'COMMON.LOCAL'
1454 include 'COMMON.NAMES'
1455 include 'COMMON.CHAIN'
1456 include 'COMMON.FFIELD'
1457 include 'COMMON.SBRIDGE'
1458 include 'COMMON.HEADER'
1459 include 'COMMON.CONTROL'
1460 include 'COMMON.DBASE'
1461 include 'COMMON.THREAD'
1462 include 'COMMON.TIME1'
1463 C Read pattern database for threading.
1464 read (icbase,*) nseq
1466 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1467 & nres_base(2,i),nres_base(3,i)
1468 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1470 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1471 c & nres_base(2,i),nres_base(3,i)
1472 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1476 if (weidis.eq.0.0D0) weidis=0.1D0
1485 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1486 write (iout,'(a,i5)') 'nexcl: ',nexcl
1487 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1490 c------------------------------------------------------------------------------
1491 subroutine setup_var
1492 implicit real*8 (a-h,o-z)
1493 include 'DIMENSIONS'
1494 include 'COMMON.IOUNITS'
1495 include 'COMMON.GEO'
1496 include 'COMMON.VAR'
1497 include 'COMMON.INTERACT'
1498 include 'COMMON.LOCAL'
1499 include 'COMMON.NAMES'
1500 include 'COMMON.CHAIN'
1501 include 'COMMON.FFIELD'
1502 include 'COMMON.SBRIDGE'
1503 include 'COMMON.HEADER'
1504 include 'COMMON.CONTROL'
1505 include 'COMMON.DBASE'
1506 include 'COMMON.THREAD'
1507 include 'COMMON.TIME1'
1508 C Set up variable list.
1514 if (itype(i).ne.10) then
1516 ialph(i,1)=nvar+nside
1520 if (indphi.gt.0) then
1522 else if (indback.gt.0) then
1527 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1530 c----------------------------------------------------------------------------
1531 subroutine gen_dist_constr
1532 C Generate CA distance constraints.
1533 implicit real*8 (a-h,o-z)
1534 include 'DIMENSIONS'
1535 include 'COMMON.IOUNITS'
1536 include 'COMMON.GEO'
1537 include 'COMMON.VAR'
1538 include 'COMMON.INTERACT'
1539 include 'COMMON.LOCAL'
1540 include 'COMMON.NAMES'
1541 include 'COMMON.CHAIN'
1542 include 'COMMON.FFIELD'
1543 include 'COMMON.SBRIDGE'
1544 include 'COMMON.HEADER'
1545 include 'COMMON.CONTROL'
1546 include 'COMMON.DBASE'
1547 include 'COMMON.THREAD'
1548 include 'COMMON.TIME1'
1549 dimension itype_pdb(maxres)
1550 common /pizda/ itype_pdb
1552 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1553 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1554 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1556 do i=nstart_sup,nstart_sup+nsup-1
1557 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1558 cd & ' seq_pdb', restyp(itype_pdb(i))
1559 do j=i+2,nstart_sup+nsup-1
1561 ihpb(nhpb)=i+nstart_seq-nstart_sup
1562 jhpb(nhpb)=j+nstart_seq-nstart_sup
1564 dhpb(nhpb)=dist(i,j)
1567 cd write (iout,'(a)') 'Distance constraints:'
1572 cd if (ii.gt.nres) then
1577 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1578 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1579 cd & dhpb(i),forcon(i)
1583 c----------------------------------------------------------------------------
1585 implicit real*8 (a-h,o-z)
1586 include 'DIMENSIONS'
1587 include 'COMMON.MAP'
1588 include 'COMMON.IOUNITS'
1589 character*3 angid(4) /'THE','PHI','ALP','OME'/
1590 character*80 mapcard,ucase
1592 read (inp,'(a)') mapcard
1593 mapcard=ucase(mapcard)
1594 if (index(mapcard,'PHI').gt.0) then
1596 else if (index(mapcard,'THE').gt.0) then
1598 else if (index(mapcard,'ALP').gt.0) then
1600 else if (index(mapcard,'OME').gt.0) then
1603 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1604 stop 'Error - illegal variable spec in MAP card.'
1606 call readi (mapcard,'RES1',res1(imap),0)
1607 call readi (mapcard,'RES2',res2(imap),0)
1608 if (res1(imap).eq.0) then
1609 res1(imap)=res2(imap)
1610 else if (res2(imap).eq.0) then
1611 res2(imap)=res1(imap)
1613 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1615 & 'Error - illegal definition of variable group in MAP.'
1616 stop 'Error - illegal definition of variable group in MAP.'
1618 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1619 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1620 call readi(mapcard,'NSTEP',nstep(imap),0)
1621 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1623 & 'Illegal boundary and/or step size specification in MAP.'
1624 stop 'Illegal boundary and/or step size specification in MAP.'
1629 c----------------------------------------------------------------------------
1630 csa subroutine csaread
1631 csa implicit real*8 (a-h,o-z)
1632 csa include 'DIMENSIONS'
1633 csa include 'COMMON.IOUNITS'
1634 csa include 'COMMON.GEO'
1635 csa include 'COMMON.CSA'
1636 csa include 'COMMON.BANK'
1637 csa include 'COMMON.CONTROL'
1638 csa character*80 ucase
1639 csa character*620 mcmcard
1640 csa call card_concat(mcmcard)
1642 csa call readi(mcmcard,'NCONF',nconf,50)
1643 csa call readi(mcmcard,'NADD',nadd,0)
1644 csa call readi(mcmcard,'JSTART',jstart,1)
1645 csa call readi(mcmcard,'JEND',jend,1)
1646 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1647 csa call readi(mcmcard,'N0',n0,1)
1648 csa call readi(mcmcard,'N1',n1,6)
1649 csa call readi(mcmcard,'N2',n2,4)
1650 csa call readi(mcmcard,'N3',n3,0)
1651 csa call readi(mcmcard,'N4',n4,0)
1652 csa call readi(mcmcard,'N5',n5,0)
1653 csa call readi(mcmcard,'N6',n6,10)
1654 csa call readi(mcmcard,'N7',n7,0)
1655 csa call readi(mcmcard,'N8',n8,0)
1656 csa call readi(mcmcard,'N9',n9,0)
1657 csa call readi(mcmcard,'N14',n14,0)
1658 csa call readi(mcmcard,'N15',n15,0)
1659 csa call readi(mcmcard,'N16',n16,0)
1660 csa call readi(mcmcard,'N17',n17,0)
1661 csa call readi(mcmcard,'N18',n18,0)
1663 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1665 csa call readi(mcmcard,'NDIFF',ndiff,2)
1666 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1667 csa call readi(mcmcard,'IS1',is1,1)
1668 csa call readi(mcmcard,'IS2',is2,8)
1669 csa call readi(mcmcard,'NRAN0',nran0,4)
1670 csa call readi(mcmcard,'NRAN1',nran1,2)
1671 csa call readi(mcmcard,'IRR',irr,1)
1672 csa call readi(mcmcard,'NSEED',nseed,20)
1673 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1674 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1675 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1676 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1677 csa call readi(mcmcard,'ICMAX',icmax,3)
1678 csa call readi(mcmcard,'IRESTART',irestart,0)
1679 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1682 csa call reada(mcmcard,'DELE',dele,20.0d0)
1683 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1684 csa call readi(mcmcard,'IREF',iref,0)
1685 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1686 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1687 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1688 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1689 csa write (iout,*) "NCONF_IN",nconf_in
1692 c----------------------------------------------------------------------------
1693 cfmc subroutine mcmfread
1694 cfmc implicit real*8 (a-h,o-z)
1695 cfmc include 'DIMENSIONS'
1696 cfmc include 'COMMON.MCMF'
1697 cfmc include 'COMMON.IOUNITS'
1698 cfmc include 'COMMON.GEO'
1699 cfmc character*80 ucase
1700 cfmc character*620 mcmcard
1701 cfmc call card_concat(mcmcard)
1703 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1704 cfmc write(iout,*)'MAXRANT=',maxrant
1705 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1706 cfmc write(iout,*)'MAXFAM=',maxfam
1707 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1708 cfmc write(iout,*)'NNET1=',nnet1
1709 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1710 cfmc write(iout,*)'NNET2=',nnet2
1711 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1712 cfmc write(iout,*)'NNET3=',nnet3
1713 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1714 cfmc write(iout,*)'ILASTT=',ilastt
1715 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1716 cfmc write(iout,*)'MAXSTR=',maxstr
1717 cfmc maxstr_f=maxstr/maxfam
1718 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1719 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1720 cfmc write(iout,*)'NMCMF=',nmcmf
1721 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1722 cfmc write(iout,*)'IFOCUS=',ifocus
1723 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1724 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1725 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1726 cfmc write(iout,*)'INTPRT=',intprt
1727 cfmc call readi(mcmcard,'IPRT',iprt,100)
1728 cfmc write(iout,*)'IPRT=',iprt
1729 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1730 cfmc write(iout,*)'IMAXTR=',imaxtr
1731 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1732 cfmc write(iout,*)'MAXEVEN=',maxeven
1733 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1734 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1735 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1736 cfmc write(iout,*)'INIMIN=',inimin
1737 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1738 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1739 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1740 cfmc write(iout,*)'NTHREAD=',nthread
1741 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1742 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1743 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1744 cfmc write(iout,*)'MAXPERT=',maxpert
1745 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1746 cfmc write(iout,*)'IRMSD=',irmsd
1747 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1748 cfmc write(iout,*)'DENEMIN=',denemin
1749 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1750 cfmc write(iout,*)'RCUT1S=',rcut1s
1751 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1752 cfmc write(iout,*)'RCUT1E=',rcut1e
1753 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1754 cfmc write(iout,*)'RCUT2S=',rcut2s
1755 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1756 cfmc write(iout,*)'RCUT2E=',rcut2e
1757 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1758 cfmc write(iout,*)'DPERT1=',d_pert1
1759 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1760 cfmc write(iout,*)'DPERT1A=',d_pert1a
1761 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1762 cfmc write(iout,*)'DPERT2=',d_pert2
1763 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1764 cfmc write(iout,*)'DPERT2A=',d_pert2a
1765 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1766 cfmc write(iout,*)'DPERT2B=',d_pert2b
1767 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1768 cfmc write(iout,*)'DPERT2C=',d_pert2c
1769 cfmc d_pert1=deg2rad*d_pert1
1770 cfmc d_pert1a=deg2rad*d_pert1a
1771 cfmc d_pert2=deg2rad*d_pert2
1772 cfmc d_pert2a=deg2rad*d_pert2a
1773 cfmc d_pert2b=deg2rad*d_pert2b
1774 cfmc d_pert2c=deg2rad*d_pert2c
1775 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1776 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1777 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1778 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1779 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1780 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1781 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1782 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1783 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1784 cfmc write(iout,*)'RCUTINI=',rcutini
1785 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1786 cfmc write(iout,*)'GRAT=',grat
1787 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1788 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1792 c----------------------------------------------------------------------------
1794 implicit real*8 (a-h,o-z)
1795 include 'DIMENSIONS'
1796 include 'COMMON.MCM'
1797 include 'COMMON.MCE'
1798 include 'COMMON.IOUNITS'
1800 character*320 mcmcard
1801 call card_concat(mcmcard)
1802 call readi(mcmcard,'MAXACC',maxacc,100)
1803 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1804 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1805 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1806 call readi(mcmcard,'MAXREPM',maxrepm,200)
1807 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1808 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1809 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1810 call reada(mcmcard,'E_UP',e_up,5.0D0)
1811 call reada(mcmcard,'DELTE',delte,0.1D0)
1812 call readi(mcmcard,'NSWEEP',nsweep,5)
1813 call readi(mcmcard,'NSTEPH',nsteph,0)
1814 call readi(mcmcard,'NSTEPC',nstepc,0)
1815 call reada(mcmcard,'TMIN',tmin,298.0D0)
1816 call reada(mcmcard,'TMAX',tmax,298.0D0)
1817 call readi(mcmcard,'NWINDOW',nwindow,0)
1818 call readi(mcmcard,'PRINT_MC',print_mc,0)
1819 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1820 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1821 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1822 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1823 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1824 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1825 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1826 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1827 if (nwindow.gt.0) then
1828 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1830 winlen(i)=winend(i)-winstart(i)+1
1833 if (tmax.lt.tmin) tmax=tmin
1834 if (tmax.eq.tmin) then
1838 if (nstepc.gt.0 .and. nsteph.gt.0) then
1839 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1840 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1842 C Probabilities of different move types
1843 sumpro_type(0)=0.0D0
1844 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1845 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1846 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1847 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1848 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1849 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1850 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1852 print *,'i',i,' sumprotype',sumpro_type(i)
1853 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1854 print *,'i',i,' sumprotype',sumpro_type(i)
1858 c----------------------------------------------------------------------------
1859 subroutine read_minim
1860 implicit real*8 (a-h,o-z)
1861 include 'DIMENSIONS'
1862 include 'COMMON.MINIM'
1863 include 'COMMON.IOUNITS'
1865 character*320 minimcard
1866 call card_concat(minimcard)
1867 call readi(minimcard,'MAXMIN',maxmin,2000)
1868 call readi(minimcard,'MAXFUN',maxfun,5000)
1869 call readi(minimcard,'MINMIN',minmin,maxmin)
1870 call readi(minimcard,'MINFUN',minfun,maxmin)
1871 call reada(minimcard,'TOLF',tolf,1.0D-2)
1872 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1873 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1874 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1875 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1876 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1877 & 'Options in energy minimization:'
1878 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1879 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1880 & 'MinMin:',MinMin,' MinFun:',MinFun,
1881 & ' TolF:',TolF,' RTolF:',RTolF
1884 c----------------------------------------------------------------------------
1885 subroutine read_angles(kanal,*)
1886 implicit real*8 (a-h,o-z)
1887 include 'DIMENSIONS'
1888 include 'COMMON.GEO'
1889 include 'COMMON.VAR'
1890 include 'COMMON.CHAIN'
1891 include 'COMMON.IOUNITS'
1892 include 'COMMON.CONTROL'
1893 c Read angles from input
1895 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1896 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1897 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1898 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1901 c 9/7/01 avoid 180 deg valence angle
1902 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1904 theta(i)=deg2rad*theta(i)
1905 phi(i)=deg2rad*phi(i)
1906 alph(i)=deg2rad*alph(i)
1907 omeg(i)=deg2rad*omeg(i)
1912 c----------------------------------------------------------------------------
1913 subroutine reada(rekord,lancuch,wartosc,default)
1915 character*(*) rekord,lancuch
1916 double precision wartosc,default
1919 iread=index(rekord,lancuch)
1920 if (iread.eq.0) then
1924 iread=iread+ilen(lancuch)+1
1925 read (rekord(iread:),*,err=10,end=10) wartosc
1930 c----------------------------------------------------------------------------
1931 subroutine readi(rekord,lancuch,wartosc,default)
1933 character*(*) rekord,lancuch
1934 integer wartosc,default
1937 iread=index(rekord,lancuch)
1938 if (iread.eq.0) then
1942 iread=iread+ilen(lancuch)+1
1943 read (rekord(iread:),*,err=10,end=10) wartosc
1948 c----------------------------------------------------------------------------
1949 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1952 integer tablica(dim),default
1953 character*(*) rekord,lancuch
1961 aux=lancuch(:ilen(lancuch))//"="
1962 iread=index(rekord,aux)
1964 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1966 if (iread.eq.0) return
1967 iread=iread+ilen(lancuch)+1
1968 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1971 c----------------------------------------------------------------------------
1972 subroutine multreada(rekord,lancuch,tablica,dim,default)
1975 double precision tablica(dim),default
1976 character*(*) rekord,lancuch
1984 aux=lancuch(:ilen(lancuch))//"="
1985 iread=index(rekord,aux)
1987 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1989 if (iread.eq.0) return
1990 iread=iread+ilen(lancuch)+1
1991 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1994 c----------------------------------------------------------------------------
1995 subroutine openunits
1996 implicit real*8 (a-h,o-z)
1997 include 'DIMENSIONS'
2000 character*16 form,nodename
2003 include 'COMMON.SETUP'
2004 include 'COMMON.IOUNITS'
2006 include 'COMMON.CONTROL'
2007 integer lenpre,lenpot,ilen,lentmp
2009 character*3 out1file_text,ucase
2012 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2013 call getenv_loc("PREFIX",prefix)
2015 call getenv_loc("POT",pot)
2016 call getenv_loc("DIRTMP",tmpdir)
2017 call getenv_loc("CURDIR",curdir)
2018 call getenv_loc("OUT1FILE",out1file_text)
2019 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2020 out1file_text=ucase(out1file_text)
2021 if (out1file_text(1:1).eq."Y") then
2024 out1file=fg_rank.gt.0
2029 if (lentmp.gt.0) then
2030 write (*,'(80(1h!))')
2031 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2032 write (*,'(80(1h!))')
2033 write (*,*)"All output files will be on node /tmp directory."
2035 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2036 if (me.eq.king) then
2037 write (*,*) "The master node is ",nodename
2038 else if (fg_rank.eq.0) then
2039 write (*,*) "I am the CG slave node ",nodename
2041 write (*,*) "I am the FG slave node ",nodename
2044 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2045 lenpre = lentmp+lenpre+1
2047 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2048 C Get the names and open the input files
2049 #if defined(WINIFL) || defined(WINPGI)
2050 open(1,file=pref_orig(:ilen(pref_orig))//
2051 & '.inp',status='old',readonly,shared)
2052 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2053 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2054 C Get parameter filenames and open the parameter files.
2055 call getenv_loc('BONDPAR',bondname)
2056 open (ibond,file=bondname,status='old',readonly,shared)
2057 call getenv_loc('THETPAR',thetname)
2058 open (ithep,file=thetname,status='old',readonly,shared)
2060 call getenv_loc('THETPARPDB',thetname_pdb)
2061 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2063 call getenv_loc('ROTPAR',rotname)
2064 open (irotam,file=rotname,status='old',readonly,shared)
2066 call getenv_loc('ROTPARPDB',rotname_pdb)
2067 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2069 call getenv_loc('TORPAR',torname)
2070 open (itorp,file=torname,status='old',readonly,shared)
2071 call getenv_loc('TORDPAR',tordname)
2072 open (itordp,file=tordname,status='old',readonly,shared)
2073 call getenv_loc('FOURIER',fouriername)
2074 open (ifourier,file=fouriername,status='old',readonly,shared)
2075 call getenv_loc('ELEPAR',elename)
2076 open (ielep,file=elename,status='old',readonly,shared)
2077 call getenv_loc('SIDEPAR',sidename)
2078 open (isidep,file=sidename,status='old',readonly,shared)
2079 #elif (defined CRAY) || (defined AIX)
2080 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2082 c print *,"Processor",myrank," opened file 1"
2083 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2084 c print *,"Processor",myrank," opened file 9"
2085 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2086 C Get parameter filenames and open the parameter files.
2087 call getenv_loc('BONDPAR',bondname)
2088 open (ibond,file=bondname,status='old',action='read')
2089 c print *,"Processor",myrank," opened file IBOND"
2090 call getenv_loc('THETPAR',thetname)
2091 open (ithep,file=thetname,status='old',action='read')
2092 c print *,"Processor",myrank," opened file ITHEP"
2094 call getenv_loc('THETPARPDB',thetname_pdb)
2095 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2097 call getenv_loc('ROTPAR',rotname)
2098 open (irotam,file=rotname,status='old',action='read')
2099 c print *,"Processor",myrank," opened file IROTAM"
2101 call getenv_loc('ROTPARPDB',rotname_pdb)
2102 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2104 call getenv_loc('TORPAR',torname)
2105 open (itorp,file=torname,status='old',action='read')
2106 c print *,"Processor",myrank," opened file ITORP"
2107 call getenv_loc('TORDPAR',tordname)
2108 open (itordp,file=tordname,status='old',action='read')
2109 c print *,"Processor",myrank," opened file ITORDP"
2110 call getenv_loc('SCCORPAR',sccorname)
2111 open (isccor,file=sccorname,status='old',action='read')
2112 c print *,"Processor",myrank," opened file ISCCOR"
2113 call getenv_loc('FOURIER',fouriername)
2114 open (ifourier,file=fouriername,status='old',action='read')
2115 c print *,"Processor",myrank," opened file IFOURIER"
2116 call getenv_loc('ELEPAR',elename)
2117 open (ielep,file=elename,status='old',action='read')
2118 c print *,"Processor",myrank," opened file IELEP"
2119 call getenv_loc('SIDEPAR',sidename)
2120 open (isidep,file=sidename,status='old',action='read')
2121 c print *,"Processor",myrank," opened file ISIDEP"
2122 c print *,"Processor",myrank," opened parameter files"
2124 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2125 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2126 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2127 C Get parameter filenames and open the parameter files.
2128 call getenv_loc('BONDPAR',bondname)
2129 open (ibond,file=bondname,status='old')
2130 call getenv_loc('THETPAR',thetname)
2131 open (ithep,file=thetname,status='old')
2133 call getenv_loc('THETPARPDB',thetname_pdb)
2134 open (ithep_pdb,file=thetname_pdb,status='old')
2136 call getenv_loc('ROTPAR',rotname)
2137 open (irotam,file=rotname,status='old')
2139 call getenv_loc('ROTPARPDB',rotname_pdb)
2140 open (irotam_pdb,file=rotname_pdb,status='old')
2142 call getenv_loc('TORPAR',torname)
2143 open (itorp,file=torname,status='old')
2144 call getenv_loc('TORDPAR',tordname)
2145 open (itordp,file=tordname,status='old')
2146 call getenv_loc('SCCORPAR',sccorname)
2147 open (isccor,file=sccorname,status='old')
2148 call getenv_loc('FOURIER',fouriername)
2149 open (ifourier,file=fouriername,status='old')
2150 call getenv_loc('ELEPAR',elename)
2151 open (ielep,file=elename,status='old')
2152 call getenv_loc('SIDEPAR',sidename)
2153 open (isidep,file=sidename,status='old')
2155 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2157 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2158 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2159 C Get parameter filenames and open the parameter files.
2160 call getenv_loc('BONDPAR',bondname)
2161 open (ibond,file=bondname,status='old',action='read')
2162 call getenv_loc('THETPAR',thetname)
2163 open (ithep,file=thetname,status='old',action='read')
2165 call getenv_loc('THETPARPDB',thetname_pdb)
2166 print *,"thetname_pdb ",thetname_pdb
2167 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2168 print *,ithep_pdb," opened"
2170 call getenv_loc('ROTPAR',rotname)
2171 open (irotam,file=rotname,status='old',action='read')
2173 call getenv_loc('ROTPARPDB',rotname_pdb)
2174 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2176 call getenv_loc('TORPAR',torname)
2177 open (itorp,file=torname,status='old',action='read')
2178 call getenv_loc('TORDPAR',tordname)
2179 open (itordp,file=tordname,status='old',action='read')
2180 call getenv_loc('SCCORPAR',sccorname)
2181 open (isccor,file=sccorname,status='old',action='read')
2182 call getenv_loc('FOURIER',fouriername)
2183 open (ifourier,file=fouriername,status='old',action='read')
2184 call getenv_loc('ELEPAR',elename)
2185 open (ielep,file=elename,status='old',action='read')
2186 call getenv_loc('SIDEPAR',sidename)
2187 open (isidep,file=sidename,status='old',action='read')
2191 C 8/9/01 In the newest version SCp interaction constants are read from a file
2192 C Use -DOLDSCP to use hard-coded constants instead.
2194 call getenv_loc('SCPPAR',scpname)
2195 #if defined(WINIFL) || defined(WINPGI)
2196 open (iscpp,file=scpname,status='old',readonly,shared)
2197 #elif (defined CRAY) || (defined AIX)
2198 open (iscpp,file=scpname,status='old',action='read')
2200 open (iscpp,file=scpname,status='old')
2202 open (iscpp,file=scpname,status='old',action='read')
2205 call getenv_loc('PATTERN',patname)
2206 #if defined(WINIFL) || defined(WINPGI)
2207 open (icbase,file=patname,status='old',readonly,shared)
2208 #elif (defined CRAY) || (defined AIX)
2209 open (icbase,file=patname,status='old',action='read')
2211 open (icbase,file=patname,status='old')
2213 open (icbase,file=patname,status='old',action='read')
2216 C Open output file only for CG processes
2217 c print *,"Processor",myrank," fg_rank",fg_rank
2218 if (fg_rank.eq.0) then
2220 if (nodes.eq.1) then
2223 npos = dlog10(dfloat(nodes-1))+1
2225 if (npos.lt.3) npos=3
2226 write (liczba,'(i1)') npos
2227 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2229 write (liczba,form) me
2230 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2231 & liczba(:ilen(liczba))
2232 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2234 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2236 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2237 & liczba(:ilen(liczba))//'.mol2'
2238 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2239 & liczba(:ilen(liczba))//'.stat'
2241 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2242 & //liczba(:ilen(liczba))//'.stat')
2243 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2246 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2247 & liczba(:ilen(liczba))//'.const'
2252 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2253 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2254 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2255 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2256 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2258 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2260 rest2name=prefix(:ilen(prefix))//'.rst'
2262 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2265 #if defined(AIX) || defined(PGI)
2266 if (me.eq.king .or. .not. out1file)
2267 & open(iout,file=outname,status='unknown')
2270 if (fg_rank.gt.0) then
2271 write (liczba,'(i3.3)') myrank/nfgtasks
2272 write (ll,'(bz,i3.3)') fg_rank
2273 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2279 open(igeom,file=intname,status='unknown',position='append')
2280 open(ipdb,file=pdbname,status='unknown')
2281 open(imol2,file=mol2name,status='unknown')
2282 open(istat,file=statname,status='unknown',position='append')
2284 c1out open(iout,file=outname,status='unknown')
2287 if (me.eq.king .or. .not.out1file)
2288 & open(iout,file=outname,status='unknown')
2291 if (fg_rank.gt.0) then
2292 print "Processor",fg_rank," opening output file"
2293 write (liczba,'(i3.3)') myrank/nfgtasks
2294 write (ll,'(bz,i3.3)') fg_rank
2295 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2301 open(igeom,file=intname,status='unknown',access='append')
2302 open(ipdb,file=pdbname,status='unknown')
2303 open(imol2,file=mol2name,status='unknown')
2304 open(istat,file=statname,status='unknown',access='append')
2306 c1out open(iout,file=outname,status='unknown')
2309 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2310 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2311 csa csa_history=prefix(:lenpre)//'.CSA.history'
2312 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2313 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2314 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2315 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2316 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2317 csa csa_int=prefix(:lenpre)//'.int'
2318 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2319 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2320 csa csa_in=prefix(:lenpre)//'.CSA.in'
2321 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2324 write (iout,'(80(1h-))')
2325 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2326 write (iout,'(80(1h-))')
2327 write (iout,*) "Input file : ",
2328 & pref_orig(:ilen(pref_orig))//'.inp'
2329 write (iout,*) "Output file : ",
2330 & outname(:ilen(outname))
2332 write (iout,*) "Sidechain potential file : ",
2333 & sidename(:ilen(sidename))
2335 write (iout,*) "SCp potential file : ",
2336 & scpname(:ilen(scpname))
2338 write (iout,*) "Electrostatic potential file : ",
2339 & elename(:ilen(elename))
2340 write (iout,*) "Cumulant coefficient file : ",
2341 & fouriername(:ilen(fouriername))
2342 write (iout,*) "Torsional parameter file : ",
2343 & torname(:ilen(torname))
2344 write (iout,*) "Double torsional parameter file : ",
2345 & tordname(:ilen(tordname))
2346 write (iout,*) "SCCOR parameter file : ",
2347 & sccorname(:ilen(sccorname))
2348 write (iout,*) "Bond & inertia constant file : ",
2349 & bondname(:ilen(bondname))
2350 write (iout,*) "Bending parameter file : ",
2351 & thetname(:ilen(thetname))
2352 write (iout,*) "Rotamer parameter file : ",
2353 & rotname(:ilen(rotname))
2354 write (iout,*) "Threading database : ",
2355 & patname(:ilen(patname))
2357 &write (iout,*)" DIRTMP : ",
2359 write (iout,'(80(1h-))')
2363 c----------------------------------------------------------------------------
2364 subroutine card_concat(card)
2365 implicit real*8 (a-h,o-z)
2366 include 'DIMENSIONS'
2367 include 'COMMON.IOUNITS'
2369 character*80 karta,ucase
2371 read (inp,'(a)') karta
2374 do while (karta(80:80).eq.'&')
2375 card=card(:ilen(card)+1)//karta(:79)
2376 read (inp,'(a)') karta
2379 card=card(:ilen(card)+1)//karta
2382 c----------------------------------------------------------------------------------
2384 implicit real*8 (a-h,o-z)
2385 include 'DIMENSIONS'
2386 include 'COMMON.CHAIN'
2387 include 'COMMON.IOUNITS'
2389 open(irest2,file=rest2name,status='unknown')
2390 read(irest2,*) totT,EK,potE,totE,t_bath
2392 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2395 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2398 read (irest2,*) iset
2403 c---------------------------------------------------------------------------------
2404 subroutine read_fragments
2405 implicit real*8 (a-h,o-z)
2406 include 'DIMENSIONS'
2410 include 'COMMON.SETUP'
2411 include 'COMMON.CHAIN'
2412 include 'COMMON.IOUNITS'
2414 include 'COMMON.CONTROL'
2415 read(inp,*) nset,nfrag,npair,nfrag_back
2416 if(me.eq.king.or..not.out1file)
2417 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2418 & " nfrag_back",nfrag_back
2420 read(inp,*) mset(iset)
2422 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2424 if(me.eq.king.or..not.out1file)
2425 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2426 & ifrag(2,i,iset), qinfrag(i,iset)
2429 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2431 if(me.eq.king.or..not.out1file)
2432 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2433 & ipair(2,i,iset), qinpair(i,iset)
2436 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2437 & wfrag_back(3,i,iset),
2438 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2439 if(me.eq.king.or..not.out1file)
2440 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2441 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2446 c-------------------------------------------------------------------------------
2447 subroutine read_dist_constr
2448 implicit real*8 (a-h,o-z)
2449 include 'DIMENSIONS'
2453 include 'COMMON.SETUP'
2454 include 'COMMON.CONTROL'
2455 include 'COMMON.CHAIN'
2456 include 'COMMON.IOUNITS'
2457 include 'COMMON.SBRIDGE'
2458 integer ifrag_(2,100),ipair_(2,100)
2459 double precision wfrag_(100),wpair_(100)
2460 character*500 controlcard
2461 c write (iout,*) "Calling read_dist_constr"
2462 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2464 call card_concat(controlcard)
2465 call readi(controlcard,"NFRAG",nfrag_,0)
2466 call readi(controlcard,"NPAIR",npair_,0)
2467 call readi(controlcard,"NDIST",ndist_,0)
2468 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2469 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2470 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2471 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2472 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2473 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2474 c write (iout,*) "IFRAG"
2476 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2478 c write (iout,*) "IPAIR"
2480 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2482 if (.not.refstr .and. nfrag.gt.0) then
2484 & "ERROR: no reference structure to compute distance restraints"
2486 & "Restraints must be specified explicitly (NDIST=number)"
2489 if (nfrag.lt.2 .and. npair.gt.0) then
2490 write (iout,*) "ERROR: Less than 2 fragments specified",
2491 & " but distance restraints between pairs requested"
2496 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2497 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2498 & ifrag_(2,i)=nstart_sup+nsup-1
2499 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2501 if (wfrag_(i).gt.0.0d0) then
2502 do j=ifrag_(1,i),ifrag_(2,i)-1
2503 do k=j+1,ifrag_(2,i)
2504 c write (iout,*) "j",j," k",k
2506 if (constr_dist.eq.1) then
2511 forcon(nhpb)=wfrag_(i)
2512 else if (constr_dist.eq.2) then
2513 if (ddjk.le.dist_cut) then
2518 forcon(nhpb)=wfrag_(i)
2525 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2528 if (.not.out1file .or. me.eq.king)
2529 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2530 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2532 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2533 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2540 if (wpair_(i).gt.0.0d0) then
2548 do j=ifrag_(1,ii),ifrag_(2,ii)
2549 do k=ifrag_(1,jj),ifrag_(2,jj)
2553 forcon(nhpb)=wpair_(i)
2554 dhpb(nhpb)=dist(j,k)
2556 if (.not.out1file .or. me.eq.king)
2557 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2558 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2560 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2561 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2568 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2569 & ibecarb(i),forcon(nhpb+1)
2570 if (forcon(nhpb+1).gt.0.0d0) then
2572 if (ibecarb(i).gt.0) then
2573 ihpb(i)=ihpb(i)+nres
2574 jhpb(i)=jhpb(i)+nres
2576 if (dhpb(nhpb).eq.0.0d0)
2577 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2581 if (.not.out1file .or. me.eq.king) then
2584 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2585 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2593 c-------------------------------------------------------------------------------
2595 subroutine flush(iu)
2600 subroutine flush(iu)
2605 c------------------------------------------------------------------------------
2606 subroutine copy_to_tmp(source)
2607 include "DIMENSIONS"
2608 include "COMMON.IOUNITS"
2609 character*(*) source
2610 character* 256 tmpfile
2614 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2615 inquire(file=tmpfile,exist=ex)
2617 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2618 & " to temporary directory..."
2619 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2620 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2624 c------------------------------------------------------------------------------
2625 subroutine move_from_tmp(source)
2626 include "DIMENSIONS"
2627 include "COMMON.IOUNITS"
2628 character*(*) source
2629 character* 256 tmpfile
2632 tmpfile=source(:ilen(source))
2633 write (*,*) "Moving ",source(:ilen(source)),
2634 & " from temporary directory to working directory"
2635 write (*,*) "/bin/mv "//tmpfile//" "//curdir
2636 call system("/bin/mv "//tmpfile//" "//curdir)
2639 c------------------------------------------------------------------------------
2640 subroutine random_init(seed)
2642 C Initialize random number generator
2644 implicit real*8 (a-h,o-z)
2645 include 'DIMENSIONS'
2651 logical OKRandom, prng_restart
2653 integer iseed_array(4)
2655 include 'COMMON.IOUNITS'
2656 include 'COMMON.TIME1'
2657 include 'COMMON.THREAD'
2658 include 'COMMON.SBRIDGE'
2659 include 'COMMON.CONTROL'
2660 include 'COMMON.MCM'
2661 include 'COMMON.MAP'
2662 include 'COMMON.HEADER'
2663 csa include 'COMMON.CSA'
2664 include 'COMMON.CHAIN'
2665 include 'COMMON.MUCA'
2667 include 'COMMON.FFIELD'
2668 include 'COMMON.SETUP'
2669 iseed=-dint(dabs(seed))
2670 if (iseed.eq.0) then
2671 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2672 & 'Random seed undefined. The program will stop.'
2673 write (*,'(/80(1h*)/20x,a/80(1h*))')
2674 & 'Random seed undefined. The program will stop.'
2676 call mpi_finalize(mpi_comm_world,ierr)
2678 stop 'Bad random seed.'
2681 if (fg_rank.eq.0) then
2685 if(me.eq.king .or. .not. out1file)
2686 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2687 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2688 OKRandom = prng_restart(me,iseedi8)
2691 tmp=65536.0d0**(4-i)
2692 iseed_array(i) = dint(seed/tmp)
2693 seed=seed-iseed_array(i)*tmp
2695 if(me.eq.king .or. .not. out1file)
2696 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2697 & (iseed_array(i),i=1,4)
2698 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2699 & (iseed_array(i),i=1,4)
2700 OKRandom = prng_restart(me,iseed_array)
2703 r1=ran_number(0.0D0,1.0D0)
2704 if(me.eq.king .or. .not. out1file)
2705 & write (iout,*) 'ran_num',r1
2706 if (r1.lt.0.0d0) OKRandom=.false.
2708 if (.not.OKRandom) then
2709 write (iout,*) 'PRNG IS NOT WORKING!!!'
2710 print *,'PRNG IS NOT WORKING!!!'
2713 call mpi_abort(mpi_comm_world,error_msg,ierr)
2716 write (iout,*) 'too many processors for parallel prng'
2717 write (*,*) 'too many processors for parallel prng'
2725 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)