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)
753 call reada(weightcard,'WDFAD',wdfa_dist,0.0d0)
754 call reada(weightcard,'WDFAT',wdfa_tor,0.0d0)
755 call reada(weightcard,'WDFAN',wdfa_nei,0.0d0)
756 call reada(weightcard,'WDFAB',wdfa_beta,0.0d0)
757 call reada(weightcard,'SCAL14',scal14,0.4D0)
758 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
759 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
760 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
761 call reada(weightcard,'TEMP0',temp0,300.0d0)
762 if (index(weightcard,'SOFT').gt.0) ipot=6
763 C 12/1/95 Added weight for the multi-body term WCORR
764 call reada(weightcard,'WCORRH',wcorr,1.0D0)
765 if (wcorr4.gt.0.0d0) wcorr=wcorr4
786 weights(24)=wdfa_dist
789 weights(27)=wdfa_beta
791 if(me.eq.king.or..not.out1file)
792 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
793 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
795 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
796 10 format (/'Energy-term weights (unscaled):'//
797 & 'WSCC= ',f10.6,' (SC-SC)'/
798 & 'WSCP= ',f10.6,' (SC-p)'/
799 & 'WELEC= ',f10.6,' (p-p electr)'/
800 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
801 & 'WBOND= ',f10.6,' (stretching)'/
802 & 'WANG= ',f10.6,' (bending)'/
803 & 'WSCLOC= ',f10.6,' (SC local)'/
804 & 'WTOR= ',f10.6,' (torsional)'/
805 & 'WTORD= ',f10.6,' (double torsional)'/
806 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
807 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
808 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
809 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
810 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
811 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
812 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
813 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
814 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
815 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
816 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
817 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
818 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
819 if(me.eq.king.or..not.out1file)then
820 if (wcorr4.gt.0.0d0) then
821 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
822 & 'between contact pairs of peptide groups'
823 write (iout,'(2(a,f5.3/))')
824 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
825 & 'Range of quenching the correlation terms:',2*delt_corr
826 else if (wcorr.gt.0.0d0) then
827 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
828 & 'between contact pairs of peptide groups'
830 write (iout,'(a,f8.3)')
831 & 'Scaling factor of 1,4 SC-p interactions:',scal14
832 write (iout,'(a,f8.3)')
833 & 'General scaling factor of SC-p interactions:',scalscp
835 r0_corr=cutoff_corr-delt_corr
837 aad(i,1)=scalscp*aad(i,1)
838 aad(i,2)=scalscp*aad(i,2)
839 bad(i,1)=scalscp*bad(i,1)
840 bad(i,2)=scalscp*bad(i,2)
842 call rescale_weights(t_bath)
843 if(me.eq.king.or..not.out1file)
844 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
845 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
847 & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta
848 22 format (/'Energy-term weights (scaled):'//
849 & 'WSCC= ',f10.6,' (SC-SC)'/
850 & 'WSCP= ',f10.6,' (SC-p)'/
851 & 'WELEC= ',f10.6,' (p-p electr)'/
852 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
853 & 'WBOND= ',f10.6,' (stretching)'/
854 & 'WANG= ',f10.6,' (bending)'/
855 & 'WSCLOC= ',f10.6,' (SC local)'/
856 & 'WTOR= ',f10.6,' (torsional)'/
857 & 'WTORD= ',f10.6,' (double torsional)'/
858 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
859 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
860 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
861 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
862 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
863 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
864 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
865 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
866 & 'WTURN6= ',f10.6,' (turns, 6th order)'/
867 & 'WDFA_D= ',f10.6,' (DFA, distance)' /
868 & 'WDFA_T= ',f10.6,' (DFA, torsional)' /
869 & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' /
870 & 'WDFA_B= ',f10.6,' (DFA, beta formation)')
872 if(me.eq.king.or..not.out1file)
873 & write (iout,*) "Reference temperature for weights calculation:",
875 call reada(weightcard,"D0CM",d0cm,3.78d0)
876 call reada(weightcard,"AKCM",akcm,15.1d0)
877 call reada(weightcard,"AKTH",akth,11.0d0)
878 call reada(weightcard,"AKCT",akct,12.0d0)
879 call reada(weightcard,"V1SS",v1ss,-1.08d0)
880 call reada(weightcard,"V2SS",v2ss,7.61d0)
881 call reada(weightcard,"V3SS",v3ss,13.7d0)
882 call reada(weightcard,"EBR",ebr,-5.50D0)
883 if(me.eq.king.or..not.out1file) then
884 write (iout,*) "Parameters of the SS-bond potential:"
885 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
887 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
888 write (iout,*) "EBR",ebr
889 print *,'indpdb=',indpdb,' pdbref=',pdbref
891 if (indpdb.gt.0 .or. pdbref) then
892 read(inp,'(a)') pdbfile
893 if(me.eq.king.or..not.out1file)
894 & write (iout,'(2a)') 'PDB data will be read from file ',
895 & pdbfile(:ilen(pdbfile))
896 open(ipdbin,file=pdbfile,status='old',err=33)
898 33 write (iout,'(a)') 'Error opening PDB file.'
901 c print *,'Begin reading pdb data'
903 c print *,'Finished reading pdb data'
904 if(me.eq.king.or..not.out1file)
905 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
906 & ' nstart_sup=',nstart_sup
908 itype_pdb(i)=itype(i)
912 nct=nstart_sup+nsup-1
913 call contact(.false.,ncont_ref,icont_ref,co)
916 C Following 2 lines for diagnostics; comment out if not needed
917 write (iout,*) "Before sideadd"
919 if(me.eq.king.or..not.out1file)
920 & write(iout,*)'Adding sidechains'
927 do while (fail.and.nsi.le.maxsi)
928 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
931 if(fail) write(iout,*)'Adding sidechain failed for res ',
932 & i,' after ',nsi,' trials'
935 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
938 C Following 2 lines for diagnostics; comment out if not needed
939 write (iout,*) "After sideadd"
942 if (indpdb.eq.0) then
943 C Read sequence if not taken from the pdb file.
945 c print *,'nres=',nres
946 if (iscode.gt.0) then
947 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
949 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
951 C Convert sequence to numeric code
953 itype(i)=rescode(i,sequence(i),iscode)
955 C Assign initial virtual bond lengths
961 vbld(i+nres)=dsc(itype(i))
962 vbld_inv(i+nres)=dsc_inv(itype(i))
963 c write (iout,*) "i",i," itype",itype(i),
964 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
968 c print '(20i4)',(itype(i),i=1,nres)
971 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
973 if (itype(i).eq.21) then
977 else if (itype(i+1).ne.20) then
979 else if (itype(i).ne.20) then
986 if(me.eq.king.or..not.out1file)then
987 write (iout,*) "ITEL"
989 write (iout,*) i,itype(i),itel(i)
991 print *,'Call Read_Bridge.'
994 C 8/13/98 Set limits to generating the dihedral angles
999 read (inp,*) ndih_constr
1000 if (ndih_constr.gt.0) then
1002 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1003 if(me.eq.king.or..not.out1file)then
1005 & 'There are',ndih_constr,' constraints on phi angles.'
1007 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1011 phi0(i)=deg2rad*phi0(i)
1012 drange(i)=deg2rad*drange(i)
1014 if(me.eq.king.or..not.out1file)
1015 & write (iout,*) 'FTORS',ftors
1018 phibound(1,ii) = phi0(i)-drange(i)
1019 phibound(2,ii) = phi0(i)+drange(i)
1024 if (me.eq.king) then
1026 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1028 write (iout,'(a3,i5,2f10.1)')
1029 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1035 cd print *,'NNT=',NNT,' NCT=',NCT
1036 if (itype(1).eq.21) nnt=2
1037 if (itype(nres).eq.21) nct=nct-1
1039 C Juyong:READ init_vars
1040 C Initialize variables!
1041 C Juyong:READ read_info
1042 C READ fragment information!!
1043 C both routines should be in dfa.F file!!
1045 if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1046 & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1048 print*, 'init_dfa_vars finished!'
1050 print*, 'read_dfa_info finished!'
1054 if(me.eq.king.or..not.out1file)
1055 & write (iout,'(a,i3)') 'nsup=',nsup
1057 if (nsup.le.(nct-nnt+1)) then
1058 do i=0,nct-nnt+1-nsup
1059 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1065 & 'Error - sequences to be superposed do not match.'
1068 do i=0,nsup-(nct-nnt+1)
1069 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1071 nstart_sup=nstart_sup+i
1077 & 'Error - sequences to be superposed do not match.'
1080 if (nsup.eq.0) nsup=nct-nnt
1081 if (nstart_sup.eq.0) nstart_sup=nnt
1082 if (nstart_seq.eq.0) nstart_seq=nnt
1083 if(me.eq.king.or..not.out1file)
1084 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1085 & ' nstart_seq=',nstart_seq
1087 c--- Zscore rms -------
1088 if (nz_start.eq.0) nz_start=nnt
1089 if (nz_end.eq.0 .and. nsup.gt.0) then
1091 else if (nz_end.eq.0) then
1094 if(me.eq.king.or..not.out1file)then
1095 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1096 write (iout,*) 'IZ_SC=',iz_sc
1098 c----------------------
1101 if (.not.pdbref) then
1102 call read_angles(inp,*38)
1104 38 write (iout,'(a)') 'Error reading reference structure.'
1106 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1107 stop 'Error reading reference structure'
1111 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1120 call contact(.true.,ncont_ref,icont_ref,co)
1122 if(me.eq.king.or..not.out1file)
1123 & write (iout,*) 'Contact order:',co
1125 if(me.eq.king.or..not.out1file)
1126 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1129 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1131 if(me.eq.king.or..not.out1file)
1132 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1133 & icont_ref(1,i),' ',
1134 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1138 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1139 if (constr_dist.gt.0) then
1140 call read_dist_constr
1142 if (nhpb.gt.0) call hpb_partition
1143 c write (iout,*) "After read_dist_constr nhpb",nhpb
1145 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1146 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1147 & modecalc.ne.10) then
1148 C If input structure hasn't been supplied from the PDB file read or generate
1150 if (iranconf.eq.0 .and. .not. extconf) then
1151 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1152 & write (iout,'(a)') 'Initial geometry will be read in.'
1154 read(inp,'(8f10.5)',end=36,err=36)
1155 & ((c(l,k),l=1,3),k=1,nres),
1156 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1157 call int_from_cart1(.false.)
1160 dc(j,i)=c(j,i+1)-c(j,i)
1161 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1165 if (itype(i).ne.10) then
1167 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1168 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1174 call read_angles(inp,*36)
1177 36 write (iout,'(a)') 'Error reading angle file.'
1179 call mpi_finalize( MPI_COMM_WORLD,IERR )
1181 stop 'Error reading angle file.'
1183 else if (extconf) then
1184 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1185 & write (iout,'(a)') 'Extended chain initial geometry.'
1187 theta(i)=90d0*deg2rad
1190 phi(i)=180d0*deg2rad
1193 alph(i)=110d0*deg2rad
1196 omeg(i)=-120d0*deg2rad
1199 if(me.eq.king.or..not.out1file)
1200 & write (iout,'(a)') 'Random-generated initial geometry.'
1204 if (me.eq.king .or. fg_rank.eq.0 .and. (
1205 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1209 call gen_rand_conf(itmp,*30)
1211 30 write (iout,*) 'Failed to generate random conformation',
1212 & ', itrial=',itrial
1213 write (*,*) 'Processor:',me,
1214 & ' Failed to generate random conformation',
1224 write (iout,'(a,i3,a)') 'Processor:',me,
1225 & ' error in generating random conformation.'
1226 write (*,'(a,i3,a)') 'Processor:',me,
1227 & ' error in generating random conformation.'
1230 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1237 elseif (modecalc.eq.4) then
1238 read (inp,'(a)') intinname
1239 open (intin,file=intinname,status='old',err=333)
1240 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1241 & write (iout,'(a)') 'intinname',intinname
1242 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1244 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1246 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1248 stop 'Error opening angle file.'
1252 C Generate distance constraints, if the PDB structure is to be regularized.
1253 if (nthread.gt.0) then
1254 call read_threadbase
1257 if (me.eq.king .or. .not. out1file)
1259 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1260 write (iout,'(/a,i3,a)')
1261 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1262 write (iout,'(20i4)') (iss(i),i=1,ns)
1263 write (iout,'(/a/)') 'Pre-formed links are:'
1269 if (me.eq.king.or..not.out1file)
1270 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1271 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1276 if (i2ndstr.gt.0) call secstrp2dihc
1277 c call geom_to_var(nvar,x)
1278 c call etotal(energia(0))
1279 c call enerprint(energia(0))
1280 c call briefout(0,etot)
1282 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1283 cd write (iout,'(a)') 'Variable list:'
1284 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1286 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1287 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1288 & 'Processor',myrank,': end reading molecular data.'
1292 c--------------------------------------------------------------------------
1293 logical function seq_comp(itypea,itypeb,length)
1295 integer length,itypea(length),itypeb(length)
1298 if (itypea(i).ne.itypeb(i)) then
1306 c-----------------------------------------------------------------------------
1307 subroutine read_bridge
1308 C Read information about disulfide bridges.
1309 implicit real*8 (a-h,o-z)
1310 include 'DIMENSIONS'
1314 include 'COMMON.IOUNITS'
1315 include 'COMMON.GEO'
1316 include 'COMMON.VAR'
1317 include 'COMMON.INTERACT'
1318 include 'COMMON.LOCAL'
1319 include 'COMMON.NAMES'
1320 include 'COMMON.CHAIN'
1321 include 'COMMON.FFIELD'
1322 include 'COMMON.SBRIDGE'
1323 include 'COMMON.HEADER'
1324 include 'COMMON.CONTROL'
1325 include 'COMMON.DBASE'
1326 include 'COMMON.THREAD'
1327 include 'COMMON.TIME1'
1328 include 'COMMON.SETUP'
1329 C Read bridging residues.
1330 read (inp,*) ns,(iss(i),i=1,ns)
1332 if(me.eq.king.or..not.out1file)
1333 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1334 C Check whether the specified bridging residues are cystines.
1336 if (itype(iss(i)).ne.1) then
1337 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1338 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1339 & ' can form a disulfide bridge?!!!'
1340 write (*,'(2a,i3,a)')
1341 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1342 & ' can form a disulfide bridge?!!!'
1344 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1349 C Read preformed bridges.
1351 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1352 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1355 C Check if the residues involved in bridges are in the specified list of
1356 C bridging residues.
1359 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1360 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1361 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1362 & ' contains residues present in other pairs.'
1363 write (*,'(a,i3,a)') 'Disulfide pair',i,
1364 & ' contains residues present in other pairs.'
1366 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1372 if (ihpb(i).eq.iss(j)) goto 10
1374 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1377 if (jhpb(i).eq.iss(j)) goto 20
1379 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1385 ihpb(i)=ihpb(i)+nres
1386 jhpb(i)=jhpb(i)+nres
1392 c----------------------------------------------------------------------------
1393 subroutine read_x(kanal,*)
1394 implicit real*8 (a-h,o-z)
1395 include 'DIMENSIONS'
1396 include 'COMMON.GEO'
1397 include 'COMMON.VAR'
1398 include 'COMMON.CHAIN'
1399 include 'COMMON.IOUNITS'
1400 include 'COMMON.CONTROL'
1401 include 'COMMON.LOCAL'
1402 include 'COMMON.INTERACT'
1403 c Read coordinates from input
1405 read(kanal,'(8f10.5)',end=10,err=10)
1406 & ((c(l,k),l=1,3),k=1,nres),
1407 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1410 c(j,2*nres)=c(j,nres)
1412 call int_from_cart1(.false.)
1415 dc(j,i)=c(j,i+1)-c(j,i)
1416 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1420 if (itype(i).ne.10) then
1422 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1423 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1431 c----------------------------------------------------------------------------
1432 subroutine read_threadbase
1433 implicit real*8 (a-h,o-z)
1434 include 'DIMENSIONS'
1435 include 'COMMON.IOUNITS'
1436 include 'COMMON.GEO'
1437 include 'COMMON.VAR'
1438 include 'COMMON.INTERACT'
1439 include 'COMMON.LOCAL'
1440 include 'COMMON.NAMES'
1441 include 'COMMON.CHAIN'
1442 include 'COMMON.FFIELD'
1443 include 'COMMON.SBRIDGE'
1444 include 'COMMON.HEADER'
1445 include 'COMMON.CONTROL'
1446 include 'COMMON.DBASE'
1447 include 'COMMON.THREAD'
1448 include 'COMMON.TIME1'
1449 C Read pattern database for threading.
1450 read (icbase,*) nseq
1452 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1453 & nres_base(2,i),nres_base(3,i)
1454 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1456 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1457 c & nres_base(2,i),nres_base(3,i)
1458 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1462 if (weidis.eq.0.0D0) weidis=0.1D0
1471 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1472 write (iout,'(a,i5)') 'nexcl: ',nexcl
1473 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1476 c------------------------------------------------------------------------------
1477 subroutine setup_var
1478 implicit real*8 (a-h,o-z)
1479 include 'DIMENSIONS'
1480 include 'COMMON.IOUNITS'
1481 include 'COMMON.GEO'
1482 include 'COMMON.VAR'
1483 include 'COMMON.INTERACT'
1484 include 'COMMON.LOCAL'
1485 include 'COMMON.NAMES'
1486 include 'COMMON.CHAIN'
1487 include 'COMMON.FFIELD'
1488 include 'COMMON.SBRIDGE'
1489 include 'COMMON.HEADER'
1490 include 'COMMON.CONTROL'
1491 include 'COMMON.DBASE'
1492 include 'COMMON.THREAD'
1493 include 'COMMON.TIME1'
1494 C Set up variable list.
1500 if (itype(i).ne.10) then
1502 ialph(i,1)=nvar+nside
1506 if (indphi.gt.0) then
1508 else if (indback.gt.0) then
1513 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1516 c----------------------------------------------------------------------------
1517 subroutine gen_dist_constr
1518 C Generate CA distance constraints.
1519 implicit real*8 (a-h,o-z)
1520 include 'DIMENSIONS'
1521 include 'COMMON.IOUNITS'
1522 include 'COMMON.GEO'
1523 include 'COMMON.VAR'
1524 include 'COMMON.INTERACT'
1525 include 'COMMON.LOCAL'
1526 include 'COMMON.NAMES'
1527 include 'COMMON.CHAIN'
1528 include 'COMMON.FFIELD'
1529 include 'COMMON.SBRIDGE'
1530 include 'COMMON.HEADER'
1531 include 'COMMON.CONTROL'
1532 include 'COMMON.DBASE'
1533 include 'COMMON.THREAD'
1534 include 'COMMON.TIME1'
1535 dimension itype_pdb(maxres)
1536 common /pizda/ itype_pdb
1538 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1539 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1540 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1542 do i=nstart_sup,nstart_sup+nsup-1
1543 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1544 cd & ' seq_pdb', restyp(itype_pdb(i))
1545 do j=i+2,nstart_sup+nsup-1
1547 ihpb(nhpb)=i+nstart_seq-nstart_sup
1548 jhpb(nhpb)=j+nstart_seq-nstart_sup
1550 dhpb(nhpb)=dist(i,j)
1553 cd write (iout,'(a)') 'Distance constraints:'
1558 cd if (ii.gt.nres) then
1563 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1564 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1565 cd & dhpb(i),forcon(i)
1569 c----------------------------------------------------------------------------
1571 implicit real*8 (a-h,o-z)
1572 include 'DIMENSIONS'
1573 include 'COMMON.MAP'
1574 include 'COMMON.IOUNITS'
1575 character*3 angid(4) /'THE','PHI','ALP','OME'/
1576 character*80 mapcard,ucase
1578 read (inp,'(a)') mapcard
1579 mapcard=ucase(mapcard)
1580 if (index(mapcard,'PHI').gt.0) then
1582 else if (index(mapcard,'THE').gt.0) then
1584 else if (index(mapcard,'ALP').gt.0) then
1586 else if (index(mapcard,'OME').gt.0) then
1589 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1590 stop 'Error - illegal variable spec in MAP card.'
1592 call readi (mapcard,'RES1',res1(imap),0)
1593 call readi (mapcard,'RES2',res2(imap),0)
1594 if (res1(imap).eq.0) then
1595 res1(imap)=res2(imap)
1596 else if (res2(imap).eq.0) then
1597 res2(imap)=res1(imap)
1599 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1601 & 'Error - illegal definition of variable group in MAP.'
1602 stop 'Error - illegal definition of variable group in MAP.'
1604 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1605 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1606 call readi(mapcard,'NSTEP',nstep(imap),0)
1607 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1609 & 'Illegal boundary and/or step size specification in MAP.'
1610 stop 'Illegal boundary and/or step size specification in MAP.'
1615 c----------------------------------------------------------------------------
1616 csa subroutine csaread
1617 csa implicit real*8 (a-h,o-z)
1618 csa include 'DIMENSIONS'
1619 csa include 'COMMON.IOUNITS'
1620 csa include 'COMMON.GEO'
1621 csa include 'COMMON.CSA'
1622 csa include 'COMMON.BANK'
1623 csa include 'COMMON.CONTROL'
1624 csa character*80 ucase
1625 csa character*620 mcmcard
1626 csa call card_concat(mcmcard)
1628 csa call readi(mcmcard,'NCONF',nconf,50)
1629 csa call readi(mcmcard,'NADD',nadd,0)
1630 csa call readi(mcmcard,'JSTART',jstart,1)
1631 csa call readi(mcmcard,'JEND',jend,1)
1632 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1633 csa call readi(mcmcard,'N0',n0,1)
1634 csa call readi(mcmcard,'N1',n1,6)
1635 csa call readi(mcmcard,'N2',n2,4)
1636 csa call readi(mcmcard,'N3',n3,0)
1637 csa call readi(mcmcard,'N4',n4,0)
1638 csa call readi(mcmcard,'N5',n5,0)
1639 csa call readi(mcmcard,'N6',n6,10)
1640 csa call readi(mcmcard,'N7',n7,0)
1641 csa call readi(mcmcard,'N8',n8,0)
1642 csa call readi(mcmcard,'N9',n9,0)
1643 csa call readi(mcmcard,'N14',n14,0)
1644 csa call readi(mcmcard,'N15',n15,0)
1645 csa call readi(mcmcard,'N16',n16,0)
1646 csa call readi(mcmcard,'N17',n17,0)
1647 csa call readi(mcmcard,'N18',n18,0)
1649 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1651 csa call readi(mcmcard,'NDIFF',ndiff,2)
1652 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1653 csa call readi(mcmcard,'IS1',is1,1)
1654 csa call readi(mcmcard,'IS2',is2,8)
1655 csa call readi(mcmcard,'NRAN0',nran0,4)
1656 csa call readi(mcmcard,'NRAN1',nran1,2)
1657 csa call readi(mcmcard,'IRR',irr,1)
1658 csa call readi(mcmcard,'NSEED',nseed,20)
1659 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1660 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1661 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1662 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1663 csa call readi(mcmcard,'ICMAX',icmax,3)
1664 csa call readi(mcmcard,'IRESTART',irestart,0)
1665 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1668 csa call reada(mcmcard,'DELE',dele,20.0d0)
1669 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1670 csa call readi(mcmcard,'IREF',iref,0)
1671 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1672 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1673 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1674 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1675 csa write (iout,*) "NCONF_IN",nconf_in
1678 c----------------------------------------------------------------------------
1679 cfmc subroutine mcmfread
1680 cfmc implicit real*8 (a-h,o-z)
1681 cfmc include 'DIMENSIONS'
1682 cfmc include 'COMMON.MCMF'
1683 cfmc include 'COMMON.IOUNITS'
1684 cfmc include 'COMMON.GEO'
1685 cfmc character*80 ucase
1686 cfmc character*620 mcmcard
1687 cfmc call card_concat(mcmcard)
1689 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1690 cfmc write(iout,*)'MAXRANT=',maxrant
1691 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1692 cfmc write(iout,*)'MAXFAM=',maxfam
1693 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1694 cfmc write(iout,*)'NNET1=',nnet1
1695 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1696 cfmc write(iout,*)'NNET2=',nnet2
1697 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1698 cfmc write(iout,*)'NNET3=',nnet3
1699 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1700 cfmc write(iout,*)'ILASTT=',ilastt
1701 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1702 cfmc write(iout,*)'MAXSTR=',maxstr
1703 cfmc maxstr_f=maxstr/maxfam
1704 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1705 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1706 cfmc write(iout,*)'NMCMF=',nmcmf
1707 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1708 cfmc write(iout,*)'IFOCUS=',ifocus
1709 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1710 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1711 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1712 cfmc write(iout,*)'INTPRT=',intprt
1713 cfmc call readi(mcmcard,'IPRT',iprt,100)
1714 cfmc write(iout,*)'IPRT=',iprt
1715 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1716 cfmc write(iout,*)'IMAXTR=',imaxtr
1717 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1718 cfmc write(iout,*)'MAXEVEN=',maxeven
1719 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1720 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1721 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1722 cfmc write(iout,*)'INIMIN=',inimin
1723 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1724 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1725 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1726 cfmc write(iout,*)'NTHREAD=',nthread
1727 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1728 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1729 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1730 cfmc write(iout,*)'MAXPERT=',maxpert
1731 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1732 cfmc write(iout,*)'IRMSD=',irmsd
1733 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1734 cfmc write(iout,*)'DENEMIN=',denemin
1735 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1736 cfmc write(iout,*)'RCUT1S=',rcut1s
1737 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1738 cfmc write(iout,*)'RCUT1E=',rcut1e
1739 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1740 cfmc write(iout,*)'RCUT2S=',rcut2s
1741 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1742 cfmc write(iout,*)'RCUT2E=',rcut2e
1743 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1744 cfmc write(iout,*)'DPERT1=',d_pert1
1745 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1746 cfmc write(iout,*)'DPERT1A=',d_pert1a
1747 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1748 cfmc write(iout,*)'DPERT2=',d_pert2
1749 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1750 cfmc write(iout,*)'DPERT2A=',d_pert2a
1751 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1752 cfmc write(iout,*)'DPERT2B=',d_pert2b
1753 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1754 cfmc write(iout,*)'DPERT2C=',d_pert2c
1755 cfmc d_pert1=deg2rad*d_pert1
1756 cfmc d_pert1a=deg2rad*d_pert1a
1757 cfmc d_pert2=deg2rad*d_pert2
1758 cfmc d_pert2a=deg2rad*d_pert2a
1759 cfmc d_pert2b=deg2rad*d_pert2b
1760 cfmc d_pert2c=deg2rad*d_pert2c
1761 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1762 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1763 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1764 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1765 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1766 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1767 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1768 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1769 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1770 cfmc write(iout,*)'RCUTINI=',rcutini
1771 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1772 cfmc write(iout,*)'GRAT=',grat
1773 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1774 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1778 c----------------------------------------------------------------------------
1780 implicit real*8 (a-h,o-z)
1781 include 'DIMENSIONS'
1782 include 'COMMON.MCM'
1783 include 'COMMON.MCE'
1784 include 'COMMON.IOUNITS'
1786 character*320 mcmcard
1787 call card_concat(mcmcard)
1788 call readi(mcmcard,'MAXACC',maxacc,100)
1789 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1790 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1791 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1792 call readi(mcmcard,'MAXREPM',maxrepm,200)
1793 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1794 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1795 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1796 call reada(mcmcard,'E_UP',e_up,5.0D0)
1797 call reada(mcmcard,'DELTE',delte,0.1D0)
1798 call readi(mcmcard,'NSWEEP',nsweep,5)
1799 call readi(mcmcard,'NSTEPH',nsteph,0)
1800 call readi(mcmcard,'NSTEPC',nstepc,0)
1801 call reada(mcmcard,'TMIN',tmin,298.0D0)
1802 call reada(mcmcard,'TMAX',tmax,298.0D0)
1803 call readi(mcmcard,'NWINDOW',nwindow,0)
1804 call readi(mcmcard,'PRINT_MC',print_mc,0)
1805 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1806 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1807 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1808 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1809 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1810 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1811 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1812 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1813 if (nwindow.gt.0) then
1814 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1816 winlen(i)=winend(i)-winstart(i)+1
1819 if (tmax.lt.tmin) tmax=tmin
1820 if (tmax.eq.tmin) then
1824 if (nstepc.gt.0 .and. nsteph.gt.0) then
1825 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1826 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1828 C Probabilities of different move types
1829 sumpro_type(0)=0.0D0
1830 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1831 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1832 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1833 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1834 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1835 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1836 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1838 print *,'i',i,' sumprotype',sumpro_type(i)
1839 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1840 print *,'i',i,' sumprotype',sumpro_type(i)
1844 c----------------------------------------------------------------------------
1845 subroutine read_minim
1846 implicit real*8 (a-h,o-z)
1847 include 'DIMENSIONS'
1848 include 'COMMON.MINIM'
1849 include 'COMMON.IOUNITS'
1851 character*320 minimcard
1852 call card_concat(minimcard)
1853 call readi(minimcard,'MAXMIN',maxmin,2000)
1854 call readi(minimcard,'MAXFUN',maxfun,5000)
1855 call readi(minimcard,'MINMIN',minmin,maxmin)
1856 call readi(minimcard,'MINFUN',minfun,maxmin)
1857 call reada(minimcard,'TOLF',tolf,1.0D-2)
1858 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1859 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1860 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1861 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1862 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1863 & 'Options in energy minimization:'
1864 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1865 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1866 & 'MinMin:',MinMin,' MinFun:',MinFun,
1867 & ' TolF:',TolF,' RTolF:',RTolF
1870 c----------------------------------------------------------------------------
1871 subroutine read_angles(kanal,*)
1872 implicit real*8 (a-h,o-z)
1873 include 'DIMENSIONS'
1874 include 'COMMON.GEO'
1875 include 'COMMON.VAR'
1876 include 'COMMON.CHAIN'
1877 include 'COMMON.IOUNITS'
1878 include 'COMMON.CONTROL'
1879 c Read angles from input
1881 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1882 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1883 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1884 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1887 c 9/7/01 avoid 180 deg valence angle
1888 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1890 theta(i)=deg2rad*theta(i)
1891 phi(i)=deg2rad*phi(i)
1892 alph(i)=deg2rad*alph(i)
1893 omeg(i)=deg2rad*omeg(i)
1898 c----------------------------------------------------------------------------
1899 subroutine reada(rekord,lancuch,wartosc,default)
1901 character*(*) rekord,lancuch
1902 double precision wartosc,default
1905 iread=index(rekord,lancuch)
1906 if (iread.eq.0) then
1910 iread=iread+ilen(lancuch)+1
1911 read (rekord(iread:),*,err=10,end=10) wartosc
1916 c----------------------------------------------------------------------------
1917 subroutine readi(rekord,lancuch,wartosc,default)
1919 character*(*) rekord,lancuch
1920 integer wartosc,default
1923 iread=index(rekord,lancuch)
1924 if (iread.eq.0) then
1928 iread=iread+ilen(lancuch)+1
1929 read (rekord(iread:),*,err=10,end=10) wartosc
1934 c----------------------------------------------------------------------------
1935 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1938 integer tablica(dim),default
1939 character*(*) rekord,lancuch
1946 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1947 if (iread.eq.0) return
1948 iread=iread+ilen(lancuch)+1
1949 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1952 c----------------------------------------------------------------------------
1953 subroutine multreada(rekord,lancuch,tablica,dim,default)
1956 double precision tablica(dim),default
1957 character*(*) rekord,lancuch
1964 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1965 if (iread.eq.0) return
1966 iread=iread+ilen(lancuch)+1
1967 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1970 c----------------------------------------------------------------------------
1971 subroutine openunits
1972 implicit real*8 (a-h,o-z)
1973 include 'DIMENSIONS'
1976 character*16 form,nodename
1979 include 'COMMON.SETUP'
1980 include 'COMMON.IOUNITS'
1982 include 'COMMON.CONTROL'
1983 integer lenpre,lenpot,ilen,lentmp
1985 character*3 out1file_text,ucase
1988 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1989 call getenv_loc("PREFIX",prefix)
1991 call getenv_loc("POT",pot)
1992 call getenv_loc("DIRTMP",tmpdir)
1993 call getenv_loc("CURDIR",curdir)
1994 call getenv_loc("OUT1FILE",out1file_text)
1995 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1996 out1file_text=ucase(out1file_text)
1997 if (out1file_text(1:1).eq."Y") then
2000 out1file=fg_rank.gt.0
2005 if (lentmp.gt.0) then
2006 write (*,'(80(1h!))')
2007 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2008 write (*,'(80(1h!))')
2009 write (*,*)"All output files will be on node /tmp directory."
2011 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2012 if (me.eq.king) then
2013 write (*,*) "The master node is ",nodename
2014 else if (fg_rank.eq.0) then
2015 write (*,*) "I am the CG slave node ",nodename
2017 write (*,*) "I am the FG slave node ",nodename
2020 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2021 lenpre = lentmp+lenpre+1
2023 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2024 C Get the names and open the input files
2025 #if defined(WINIFL) || defined(WINPGI)
2026 open(1,file=pref_orig(:ilen(pref_orig))//
2027 & '.inp',status='old',readonly,shared)
2028 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2029 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2030 C Get parameter filenames and open the parameter files.
2031 call getenv_loc('BONDPAR',bondname)
2032 open (ibond,file=bondname,status='old',readonly,shared)
2033 call getenv_loc('THETPAR',thetname)
2034 open (ithep,file=thetname,status='old',readonly,shared)
2036 call getenv_loc('THETPARPDB',thetname_pdb)
2037 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2039 call getenv_loc('ROTPAR',rotname)
2040 open (irotam,file=rotname,status='old',readonly,shared)
2042 call getenv_loc('ROTPARPDB',rotname_pdb)
2043 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2045 call getenv_loc('TORPAR',torname)
2046 open (itorp,file=torname,status='old',readonly,shared)
2047 call getenv_loc('TORDPAR',tordname)
2048 open (itordp,file=tordname,status='old',readonly,shared)
2049 call getenv_loc('FOURIER',fouriername)
2050 open (ifourier,file=fouriername,status='old',readonly,shared)
2051 call getenv_loc('ELEPAR',elename)
2052 open (ielep,file=elename,status='old',readonly,shared)
2053 call getenv_loc('SIDEPAR',sidename)
2054 open (isidep,file=sidename,status='old',readonly,shared)
2055 #elif (defined CRAY) || (defined AIX)
2056 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2058 c print *,"Processor",myrank," opened file 1"
2059 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2060 c print *,"Processor",myrank," opened file 9"
2061 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2062 C Get parameter filenames and open the parameter files.
2063 call getenv_loc('BONDPAR',bondname)
2064 open (ibond,file=bondname,status='old',action='read')
2065 c print *,"Processor",myrank," opened file IBOND"
2066 call getenv_loc('THETPAR',thetname)
2067 open (ithep,file=thetname,status='old',action='read')
2068 c print *,"Processor",myrank," opened file ITHEP"
2070 call getenv_loc('THETPARPDB',thetname_pdb)
2071 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2073 call getenv_loc('ROTPAR',rotname)
2074 open (irotam,file=rotname,status='old',action='read')
2075 c print *,"Processor",myrank," opened file IROTAM"
2077 call getenv_loc('ROTPARPDB',rotname_pdb)
2078 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2080 call getenv_loc('TORPAR',torname)
2081 open (itorp,file=torname,status='old',action='read')
2082 c print *,"Processor",myrank," opened file ITORP"
2083 call getenv_loc('TORDPAR',tordname)
2084 open (itordp,file=tordname,status='old',action='read')
2085 c print *,"Processor",myrank," opened file ITORDP"
2086 call getenv_loc('SCCORPAR',sccorname)
2087 open (isccor,file=sccorname,status='old',action='read')
2088 c print *,"Processor",myrank," opened file ISCCOR"
2089 call getenv_loc('FOURIER',fouriername)
2090 open (ifourier,file=fouriername,status='old',action='read')
2091 c print *,"Processor",myrank," opened file IFOURIER"
2092 call getenv_loc('ELEPAR',elename)
2093 open (ielep,file=elename,status='old',action='read')
2094 c print *,"Processor",myrank," opened file IELEP"
2095 call getenv_loc('SIDEPAR',sidename)
2096 open (isidep,file=sidename,status='old',action='read')
2097 c print *,"Processor",myrank," opened file ISIDEP"
2098 c print *,"Processor",myrank," opened parameter files"
2100 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2101 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2102 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2103 C Get parameter filenames and open the parameter files.
2104 call getenv_loc('BONDPAR',bondname)
2105 open (ibond,file=bondname,status='old')
2106 call getenv_loc('THETPAR',thetname)
2107 open (ithep,file=thetname,status='old')
2109 call getenv_loc('THETPARPDB',thetname_pdb)
2110 open (ithep_pdb,file=thetname_pdb,status='old')
2112 call getenv_loc('ROTPAR',rotname)
2113 open (irotam,file=rotname,status='old')
2115 call getenv_loc('ROTPARPDB',rotname_pdb)
2116 open (irotam_pdb,file=rotname_pdb,status='old')
2118 call getenv_loc('TORPAR',torname)
2119 open (itorp,file=torname,status='old')
2120 call getenv_loc('TORDPAR',tordname)
2121 open (itordp,file=tordname,status='old')
2122 call getenv_loc('SCCORPAR',sccorname)
2123 open (isccor,file=sccorname,status='old')
2124 call getenv_loc('FOURIER',fouriername)
2125 open (ifourier,file=fouriername,status='old')
2126 call getenv_loc('ELEPAR',elename)
2127 open (ielep,file=elename,status='old')
2128 call getenv_loc('SIDEPAR',sidename)
2129 open (isidep,file=sidename,status='old')
2131 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2133 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2134 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2135 C Get parameter filenames and open the parameter files.
2136 call getenv_loc('BONDPAR',bondname)
2137 open (ibond,file=bondname,status='old',action='read')
2138 call getenv_loc('THETPAR',thetname)
2139 open (ithep,file=thetname,status='old',action='read')
2141 call getenv_loc('THETPARPDB',thetname_pdb)
2142 print *,"thetname_pdb ",thetname_pdb
2143 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2144 print *,ithep_pdb," opened"
2146 call getenv_loc('ROTPAR',rotname)
2147 open (irotam,file=rotname,status='old',action='read')
2149 call getenv_loc('ROTPARPDB',rotname_pdb)
2150 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2152 call getenv_loc('TORPAR',torname)
2153 open (itorp,file=torname,status='old',action='read')
2154 call getenv_loc('TORDPAR',tordname)
2155 open (itordp,file=tordname,status='old',action='read')
2156 call getenv_loc('SCCORPAR',sccorname)
2157 open (isccor,file=sccorname,status='old',action='read')
2158 call getenv_loc('FOURIER',fouriername)
2159 open (ifourier,file=fouriername,status='old',action='read')
2160 call getenv_loc('ELEPAR',elename)
2161 open (ielep,file=elename,status='old',action='read')
2162 call getenv_loc('SIDEPAR',sidename)
2163 open (isidep,file=sidename,status='old',action='read')
2167 C 8/9/01 In the newest version SCp interaction constants are read from a file
2168 C Use -DOLDSCP to use hard-coded constants instead.
2170 call getenv_loc('SCPPAR',scpname)
2171 #if defined(WINIFL) || defined(WINPGI)
2172 open (iscpp,file=scpname,status='old',readonly,shared)
2173 #elif (defined CRAY) || (defined AIX)
2174 open (iscpp,file=scpname,status='old',action='read')
2176 open (iscpp,file=scpname,status='old')
2178 open (iscpp,file=scpname,status='old',action='read')
2181 call getenv_loc('PATTERN',patname)
2182 #if defined(WINIFL) || defined(WINPGI)
2183 open (icbase,file=patname,status='old',readonly,shared)
2184 #elif (defined CRAY) || (defined AIX)
2185 open (icbase,file=patname,status='old',action='read')
2187 open (icbase,file=patname,status='old')
2189 open (icbase,file=patname,status='old',action='read')
2192 C Open output file only for CG processes
2193 c print *,"Processor",myrank," fg_rank",fg_rank
2194 if (fg_rank.eq.0) then
2196 if (nodes.eq.1) then
2199 npos = dlog10(dfloat(nodes-1))+1
2201 if (npos.lt.3) npos=3
2202 write (liczba,'(i1)') npos
2203 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2205 write (liczba,form) me
2206 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2207 & liczba(:ilen(liczba))
2208 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2210 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2212 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2213 & liczba(:ilen(liczba))//'.mol2'
2214 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2215 & liczba(:ilen(liczba))//'.stat'
2217 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2218 & //liczba(:ilen(liczba))//'.stat')
2219 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2222 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2223 & liczba(:ilen(liczba))//'.const'
2228 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2229 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2230 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2231 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2232 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2234 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2236 rest2name=prefix(:ilen(prefix))//'.rst'
2238 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2241 #if defined(AIX) || defined(PGI)
2242 if (me.eq.king .or. .not. out1file)
2243 & open(iout,file=outname,status='unknown')
2246 if (fg_rank.gt.0) then
2247 write (liczba,'(i3.3)') myrank/nfgtasks
2248 write (ll,'(bz,i3.3)') fg_rank
2249 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2255 open(igeom,file=intname,status='unknown',position='append')
2256 open(ipdb,file=pdbname,status='unknown')
2257 open(imol2,file=mol2name,status='unknown')
2258 open(istat,file=statname,status='unknown',position='append')
2260 c1out open(iout,file=outname,status='unknown')
2263 if (me.eq.king .or. .not.out1file)
2264 & open(iout,file=outname,status='unknown')
2267 if (fg_rank.gt.0) then
2268 print "Processor",fg_rank," opening output file"
2269 write (liczba,'(i3.3)') myrank/nfgtasks
2270 write (ll,'(bz,i3.3)') fg_rank
2271 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2277 open(igeom,file=intname,status='unknown',access='append')
2278 open(ipdb,file=pdbname,status='unknown')
2279 open(imol2,file=mol2name,status='unknown')
2280 open(istat,file=statname,status='unknown',access='append')
2282 c1out open(iout,file=outname,status='unknown')
2285 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2286 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2287 csa csa_history=prefix(:lenpre)//'.CSA.history'
2288 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2289 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2290 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2291 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2292 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2293 csa csa_int=prefix(:lenpre)//'.int'
2294 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2295 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2296 csa csa_in=prefix(:lenpre)//'.CSA.in'
2297 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2300 write (iout,'(80(1h-))')
2301 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2302 write (iout,'(80(1h-))')
2303 write (iout,*) "Input file : ",
2304 & pref_orig(:ilen(pref_orig))//'.inp'
2305 write (iout,*) "Output file : ",
2306 & outname(:ilen(outname))
2308 write (iout,*) "Sidechain potential file : ",
2309 & sidename(:ilen(sidename))
2311 write (iout,*) "SCp potential file : ",
2312 & scpname(:ilen(scpname))
2314 write (iout,*) "Electrostatic potential file : ",
2315 & elename(:ilen(elename))
2316 write (iout,*) "Cumulant coefficient file : ",
2317 & fouriername(:ilen(fouriername))
2318 write (iout,*) "Torsional parameter file : ",
2319 & torname(:ilen(torname))
2320 write (iout,*) "Double torsional parameter file : ",
2321 & tordname(:ilen(tordname))
2322 write (iout,*) "SCCOR parameter file : ",
2323 & sccorname(:ilen(sccorname))
2324 write (iout,*) "Bond & inertia constant file : ",
2325 & bondname(:ilen(bondname))
2326 write (iout,*) "Bending parameter file : ",
2327 & thetname(:ilen(thetname))
2328 write (iout,*) "Rotamer parameter file : ",
2329 & rotname(:ilen(rotname))
2330 write (iout,*) "Threading database : ",
2331 & patname(:ilen(patname))
2333 &write (iout,*)" DIRTMP : ",
2335 write (iout,'(80(1h-))')
2339 c----------------------------------------------------------------------------
2340 subroutine card_concat(card)
2341 implicit real*8 (a-h,o-z)
2342 include 'DIMENSIONS'
2343 include 'COMMON.IOUNITS'
2345 character*80 karta,ucase
2347 read (inp,'(a)') karta
2350 do while (karta(80:80).eq.'&')
2351 card=card(:ilen(card)+1)//karta(:79)
2352 read (inp,'(a)') karta
2355 card=card(:ilen(card)+1)//karta
2358 c----------------------------------------------------------------------------------
2360 implicit real*8 (a-h,o-z)
2361 include 'DIMENSIONS'
2362 include 'COMMON.CHAIN'
2363 include 'COMMON.IOUNITS'
2365 open(irest2,file=rest2name,status='unknown')
2366 read(irest2,*) totT,EK,potE,totE,t_bath
2368 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2371 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2374 read (irest2,*) iset
2379 c---------------------------------------------------------------------------------
2380 subroutine read_fragments
2381 implicit real*8 (a-h,o-z)
2382 include 'DIMENSIONS'
2386 include 'COMMON.SETUP'
2387 include 'COMMON.CHAIN'
2388 include 'COMMON.IOUNITS'
2390 include 'COMMON.CONTROL'
2391 read(inp,*) nset,nfrag,npair,nfrag_back
2392 if(me.eq.king.or..not.out1file)
2393 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2394 & " nfrag_back",nfrag_back
2396 read(inp,*) mset(iset)
2398 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2400 if(me.eq.king.or..not.out1file)
2401 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2402 & ifrag(2,i,iset), qinfrag(i,iset)
2405 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2407 if(me.eq.king.or..not.out1file)
2408 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2409 & ipair(2,i,iset), qinpair(i,iset)
2412 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2413 & wfrag_back(3,i,iset),
2414 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2415 if(me.eq.king.or..not.out1file)
2416 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2417 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2422 c-------------------------------------------------------------------------------
2423 subroutine read_dist_constr
2424 implicit real*8 (a-h,o-z)
2425 include 'DIMENSIONS'
2429 include 'COMMON.SETUP'
2430 include 'COMMON.CONTROL'
2431 include 'COMMON.CHAIN'
2432 include 'COMMON.IOUNITS'
2433 include 'COMMON.SBRIDGE'
2434 integer ifrag_(2,100),ipair_(2,100)
2435 double precision wfrag_(100),wpair_(100)
2436 character*500 controlcard
2437 c write (iout,*) "Calling read_dist_constr"
2438 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2440 call card_concat(controlcard)
2441 call readi(controlcard,"NFRAG",nfrag_,0)
2442 call readi(controlcard,"NPAIR",npair_,0)
2443 call readi(controlcard,"NDIST",ndist_,0)
2444 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2445 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2446 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2447 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2448 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2449 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2450 c write (iout,*) "IFRAG"
2452 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2454 c write (iout,*) "IPAIR"
2456 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2458 if (.not.refstr .and. nfrag.gt.0) then
2460 & "ERROR: no reference structure to compute distance restraints"
2462 & "Restraints must be specified explicitly (NDIST=number)"
2465 if (nfrag.lt.2 .and. npair.gt.0) then
2466 write (iout,*) "ERROR: Less than 2 fragments specified",
2467 & " but distance restraints between pairs requested"
2472 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2473 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2474 & ifrag_(2,i)=nstart_sup+nsup-1
2475 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2477 if (wfrag_(i).gt.0.0d0) then
2478 do j=ifrag_(1,i),ifrag_(2,i)-1
2479 do k=j+1,ifrag_(2,i)
2480 write (iout,*) "j",j," k",k
2482 if (constr_dist.eq.1) then
2487 forcon(nhpb)=wfrag_(i)
2488 else if (constr_dist.eq.2) then
2489 if (ddjk.le.dist_cut) then
2494 forcon(nhpb)=wfrag_(i)
2501 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2504 if (.not.out1file .or. me.eq.king)
2505 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2506 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2508 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2509 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2516 if (wpair_(i).gt.0.0d0) then
2524 do j=ifrag_(1,ii),ifrag_(2,ii)
2525 do k=ifrag_(1,jj),ifrag_(2,jj)
2529 forcon(nhpb)=wpair_(i)
2530 dhpb(nhpb)=dist(j,k)
2532 if (.not.out1file .or. me.eq.king)
2533 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2534 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2536 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2537 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2544 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2545 & ibecarb(i),forcon(nhpb+1)
2546 if (forcon(nhpb+1).gt.0.0d0) then
2548 if (ibecarb(i).gt.0) then
2549 ihpb(i)=ihpb(i)+nres
2550 jhpb(i)=jhpb(i)+nres
2552 if (dhpb(nhpb).eq.0.0d0)
2553 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2557 if (.not.out1file .or. me.eq.king) then
2560 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2561 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2569 c-------------------------------------------------------------------------------
2571 subroutine flush(iu)
2576 subroutine flush(iu)
2581 c------------------------------------------------------------------------------
2582 subroutine copy_to_tmp(source)
2583 include "DIMENSIONS"
2584 include "COMMON.IOUNITS"
2585 character*(*) source
2586 character* 256 tmpfile
2590 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2591 inquire(file=tmpfile,exist=ex)
2593 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2594 & " to temporary directory..."
2595 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2596 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2600 c------------------------------------------------------------------------------
2601 subroutine move_from_tmp(source)
2602 include "DIMENSIONS"
2603 include "COMMON.IOUNITS"
2604 character*(*) source
2607 write (*,*) "Moving ",source(:ilen(source)),
2608 & " from temporary directory to working directory"
2609 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2610 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2613 c------------------------------------------------------------------------------
2614 subroutine random_init(seed)
2616 C Initialize random number generator
2618 implicit real*8 (a-h,o-z)
2619 include 'DIMENSIONS'
2625 logical OKRandom, prng_restart
2627 integer iseed_array(4)
2629 include 'COMMON.IOUNITS'
2630 include 'COMMON.TIME1'
2631 include 'COMMON.THREAD'
2632 include 'COMMON.SBRIDGE'
2633 include 'COMMON.CONTROL'
2634 include 'COMMON.MCM'
2635 include 'COMMON.MAP'
2636 include 'COMMON.HEADER'
2637 csa include 'COMMON.CSA'
2638 include 'COMMON.CHAIN'
2639 include 'COMMON.MUCA'
2641 include 'COMMON.FFIELD'
2642 include 'COMMON.SETUP'
2643 iseed=-dint(dabs(seed))
2644 if (iseed.eq.0) then
2645 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2646 & 'Random seed undefined. The program will stop.'
2647 write (*,'(/80(1h*)/20x,a/80(1h*))')
2648 & 'Random seed undefined. The program will stop.'
2650 call mpi_finalize(mpi_comm_world,ierr)
2652 stop 'Bad random seed.'
2655 if (fg_rank.eq.0) then
2659 if(me.eq.king .or. .not. out1file)
2660 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2661 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2662 OKRandom = prng_restart(me,iseedi8)
2665 tmp=65536.0d0**(4-i)
2666 iseed_array(i) = dint(seed/tmp)
2667 seed=seed-iseed_array(i)*tmp
2669 if(me.eq.king .or. .not. out1file)
2670 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2671 & (iseed_array(i),i=1,4)
2672 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2673 & (iseed_array(i),i=1,4)
2674 OKRandom = prng_restart(me,iseed_array)
2677 r1=ran_number(0.0D0,1.0D0)
2678 if(me.eq.king .or. .not. out1file)
2679 & write (iout,*) 'ran_num',r1
2680 if (r1.lt.0.0d0) OKRandom=.false.
2682 if (.not.OKRandom) then
2683 write (iout,*) 'PRNG IS NOT WORKING!!!'
2684 print *,'PRNG IS NOT WORKING!!!'
2687 call mpi_abort(mpi_comm_world,error_msg,ierr)
2690 write (iout,*) 'too many processors for parallel prng'
2691 write (*,*) 'too many processors for parallel prng'
2699 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)