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 if(me.eq.king.or..not.out1file) then
864 write (iout,*) "Parameters of the SS-bond potential:"
865 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
867 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
868 write (iout,*) "EBR",ebr
869 print *,'indpdb=',indpdb,' pdbref=',pdbref
871 if (indpdb.gt.0 .or. pdbref) then
872 read(inp,'(a)') pdbfile
873 if(me.eq.king.or..not.out1file)
874 & write (iout,'(2a)') 'PDB data will be read from file ',
875 & pdbfile(:ilen(pdbfile))
876 open(ipdbin,file=pdbfile,status='old',err=33)
878 33 write (iout,'(a)') 'Error opening PDB file.'
881 c print *,'Begin reading pdb data'
883 c print *,'Finished reading pdb data'
884 if(me.eq.king.or..not.out1file)
885 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
886 & ' nstart_sup=',nstart_sup
888 itype_pdb(i)=itype(i)
892 nct=nstart_sup+nsup-1
893 call contact(.false.,ncont_ref,icont_ref,co)
896 if(me.eq.king.or..not.out1file)
897 & write(iout,*)'Adding sidechains'
904 do while (fail.and.nsi.le.maxsi)
905 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
908 if(fail) write(iout,*)'Adding sidechain failed for res ',
909 & i,' after ',nsi,' trials'
914 if (indpdb.eq.0) then
915 C Read sequence if not taken from the pdb file.
917 c print *,'nres=',nres
918 if (iscode.gt.0) then
919 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
921 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
923 C Convert sequence to numeric code
925 itype(i)=rescode(i,sequence(i),iscode)
927 C Assign initial virtual bond lengths
933 vbld(i+nres)=dsc(iabs(itype(i)))
934 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
935 c write (iout,*) "i",i," itype",itype(i),
936 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
940 c print '(20i4)',(itype(i),i=1,nres)
943 if (itype(i).eq.21 .or. itype(i+1).eq.21) then
945 if (itype(i).eq.21) then
949 else if (iabs(itype(i+1)).ne.20) then
951 else if (iabs(itype(i)).ne.20) then
958 if(me.eq.king.or..not.out1file)then
959 write (iout,*) "ITEL"
961 write (iout,*) i,itype(i),itel(i)
963 print *,'Call Read_Bridge.'
966 C 8/13/98 Set limits to generating the dihedral angles
971 read (inp,*) ndih_constr
972 if (ndih_constr.gt.0) then
974 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
975 if(me.eq.king.or..not.out1file)then
977 & 'There are',ndih_constr,' constraints on phi angles.'
979 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
983 phi0(i)=deg2rad*phi0(i)
984 drange(i)=deg2rad*drange(i)
986 if(me.eq.king.or..not.out1file)
987 & write (iout,*) 'FTORS',ftors
990 phibound(1,ii) = phi0(i)-drange(i)
991 phibound(2,ii) = phi0(i)+drange(i)
998 write (iout,'(a)') 'Boundaries in phi angle sampling:'
1000 write (iout,'(a3,i5,2f10.1)')
1001 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1007 cd print *,'NNT=',NNT,' NCT=',NCT
1008 if (itype(1).eq.21) nnt=2
1009 if (itype(nres).eq.21) nct=nct-1
1011 if(me.eq.king.or..not.out1file)
1012 & write (iout,'(a,i3)') 'nsup=',nsup
1014 if (nsup.le.(nct-nnt+1)) then
1015 do i=0,nct-nnt+1-nsup
1016 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1022 & 'Error - sequences to be superposed do not match.'
1025 do i=0,nsup-(nct-nnt+1)
1026 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1028 nstart_sup=nstart_sup+i
1034 & 'Error - sequences to be superposed do not match.'
1037 if (nsup.eq.0) nsup=nct-nnt
1038 if (nstart_sup.eq.0) nstart_sup=nnt
1039 if (nstart_seq.eq.0) nstart_seq=nnt
1040 if(me.eq.king.or..not.out1file)
1041 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1042 & ' nstart_seq=',nstart_seq
1044 c--- Zscore rms -------
1045 if (nz_start.eq.0) nz_start=nnt
1046 if (nz_end.eq.0 .and. nsup.gt.0) then
1048 else if (nz_end.eq.0) then
1051 if(me.eq.king.or..not.out1file)then
1052 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1053 write (iout,*) 'IZ_SC=',iz_sc
1055 c----------------------
1058 if (.not.pdbref) then
1059 call read_angles(inp,*38)
1061 38 write (iout,'(a)') 'Error reading reference structure.'
1063 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1064 stop 'Error reading reference structure'
1068 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
1077 call contact(.true.,ncont_ref,icont_ref,co)
1079 if(me.eq.king.or..not.out1file)
1080 & write (iout,*) 'Contact order:',co
1082 if(me.eq.king.or..not.out1file)
1083 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1086 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1088 if(me.eq.king.or..not.out1file)
1089 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1090 & icont_ref(1,i),' ',
1091 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1095 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1096 if (constr_dist.gt.0) then
1097 call read_dist_constr
1100 c write (iout,*) "After read_dist_constr nhpb",nhpb
1102 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1103 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1104 & modecalc.ne.10) then
1105 C If input structure hasn't been supplied from the PDB file read or generate
1107 if (iranconf.eq.0 .and. .not. extconf) then
1108 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1109 & write (iout,'(a)') 'Initial geometry will be read in.'
1111 read(inp,'(8f10.5)',end=36,err=36)
1112 & ((c(l,k),l=1,3),k=1,nres),
1113 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1114 call int_from_cart1(.false.)
1117 dc(j,i)=c(j,i+1)-c(j,i)
1118 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1122 if (itype(i).ne.10) then
1124 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1125 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1131 call read_angles(inp,*36)
1134 36 write (iout,'(a)') 'Error reading angle file.'
1136 call mpi_finalize( MPI_COMM_WORLD,IERR )
1138 stop 'Error reading angle file.'
1140 else if (extconf) then
1141 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1142 & write (iout,'(a)') 'Extended chain initial geometry.'
1144 theta(i)=90d0*deg2rad
1147 phi(i)=180d0*deg2rad
1150 alph(i)=110d0*deg2rad
1153 omeg(i)=-120d0*deg2rad
1154 if (itype(i).le.0) omeg(i)=-omeg(i)
1157 if(me.eq.king.or..not.out1file)
1158 & write (iout,'(a)') 'Random-generated initial geometry.'
1162 if (me.eq.king .or. fg_rank.eq.0 .and. (
1163 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1167 call gen_rand_conf(itmp,*30)
1169 30 write (iout,*) 'Failed to generate random conformation',
1170 & ', itrial=',itrial
1171 write (*,*) 'Processor:',me,
1172 & ' Failed to generate random conformation',
1182 write (iout,'(a,i3,a)') 'Processor:',me,
1183 & ' error in generating random conformation.'
1184 write (*,'(a,i3,a)') 'Processor:',me,
1185 & ' error in generating random conformation.'
1188 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1195 elseif (modecalc.eq.4) then
1196 read (inp,'(a)') intinname
1197 open (intin,file=intinname,status='old',err=333)
1198 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1199 & write (iout,'(a)') 'intinname',intinname
1200 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1202 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1204 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1206 stop 'Error opening angle file.'
1210 C Generate distance constraints, if the PDB structure is to be regularized.
1211 if (nthread.gt.0) then
1212 call read_threadbase
1215 if (me.eq.king .or. .not. out1file)
1217 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1218 write (iout,'(/a,i3,a)')
1219 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1220 write (iout,'(20i4)') (iss(i),i=1,ns)
1221 write (iout,'(/a/)') 'Pre-formed links are:'
1227 if (me.eq.king.or..not.out1file)
1228 & write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1229 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1234 if (i2ndstr.gt.0) call secstrp2dihc
1235 c call geom_to_var(nvar,x)
1236 c call etotal(energia(0))
1237 c call enerprint(energia(0))
1238 c call briefout(0,etot)
1240 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1241 cd write (iout,'(a)') 'Variable list:'
1242 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1244 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1245 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1246 & 'Processor',myrank,': end reading molecular data.'
1250 c--------------------------------------------------------------------------
1251 logical function seq_comp(itypea,itypeb,length)
1253 integer length,itypea(length),itypeb(length)
1256 if (itypea(i).ne.itypeb(i)) then
1264 c-----------------------------------------------------------------------------
1265 subroutine read_bridge
1266 C Read information about disulfide bridges.
1267 implicit real*8 (a-h,o-z)
1268 include 'DIMENSIONS'
1272 include 'COMMON.IOUNITS'
1273 include 'COMMON.GEO'
1274 include 'COMMON.VAR'
1275 include 'COMMON.INTERACT'
1276 include 'COMMON.LOCAL'
1277 include 'COMMON.NAMES'
1278 include 'COMMON.CHAIN'
1279 include 'COMMON.FFIELD'
1280 include 'COMMON.SBRIDGE'
1281 include 'COMMON.HEADER'
1282 include 'COMMON.CONTROL'
1283 include 'COMMON.DBASE'
1284 include 'COMMON.THREAD'
1285 include 'COMMON.TIME1'
1286 include 'COMMON.SETUP'
1287 C Read bridging residues.
1288 read (inp,*) ns,(iss(i),i=1,ns)
1290 if(me.eq.king.or..not.out1file)
1291 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1292 C Check whether the specified bridging residues are cystines.
1294 if (itype(iss(i)).ne.1) then
1295 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1296 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1297 & ' can form a disulfide bridge?!!!'
1298 write (*,'(2a,i3,a)')
1299 & 'Do you REALLY think that the residue ',restyp(iss(i)),i,
1300 & ' can form a disulfide bridge?!!!'
1302 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1307 C Read preformed bridges.
1309 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1310 write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1313 C Check if the residues involved in bridges are in the specified list of
1314 C bridging residues.
1317 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1318 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1319 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1320 & ' contains residues present in other pairs.'
1321 write (*,'(a,i3,a)') 'Disulfide pair',i,
1322 & ' contains residues present in other pairs.'
1324 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1330 if (ihpb(i).eq.iss(j)) goto 10
1332 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1335 if (jhpb(i).eq.iss(j)) goto 20
1337 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1343 ihpb(i)=ihpb(i)+nres
1344 jhpb(i)=jhpb(i)+nres
1350 c----------------------------------------------------------------------------
1351 subroutine read_x(kanal,*)
1352 implicit real*8 (a-h,o-z)
1353 include 'DIMENSIONS'
1354 include 'COMMON.GEO'
1355 include 'COMMON.VAR'
1356 include 'COMMON.CHAIN'
1357 include 'COMMON.IOUNITS'
1358 include 'COMMON.CONTROL'
1359 include 'COMMON.LOCAL'
1360 include 'COMMON.INTERACT'
1361 c Read coordinates from input
1363 read(kanal,'(8f10.5)',end=10,err=10)
1364 & ((c(l,k),l=1,3),k=1,nres),
1365 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1368 c(j,2*nres)=c(j,nres)
1370 call int_from_cart1(.false.)
1373 dc(j,i)=c(j,i+1)-c(j,i)
1374 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1378 if (itype(i).ne.10) then
1380 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1381 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1389 c----------------------------------------------------------------------------
1390 subroutine read_threadbase
1391 implicit real*8 (a-h,o-z)
1392 include 'DIMENSIONS'
1393 include 'COMMON.IOUNITS'
1394 include 'COMMON.GEO'
1395 include 'COMMON.VAR'
1396 include 'COMMON.INTERACT'
1397 include 'COMMON.LOCAL'
1398 include 'COMMON.NAMES'
1399 include 'COMMON.CHAIN'
1400 include 'COMMON.FFIELD'
1401 include 'COMMON.SBRIDGE'
1402 include 'COMMON.HEADER'
1403 include 'COMMON.CONTROL'
1404 include 'COMMON.DBASE'
1405 include 'COMMON.THREAD'
1406 include 'COMMON.TIME1'
1407 C Read pattern database for threading.
1408 read (icbase,*) nseq
1410 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1411 & nres_base(2,i),nres_base(3,i)
1412 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1414 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1415 c & nres_base(2,i),nres_base(3,i)
1416 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1420 if (weidis.eq.0.0D0) weidis=0.1D0
1429 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1430 write (iout,'(a,i5)') 'nexcl: ',nexcl
1431 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1434 c------------------------------------------------------------------------------
1435 subroutine setup_var
1436 implicit real*8 (a-h,o-z)
1437 include 'DIMENSIONS'
1438 include 'COMMON.IOUNITS'
1439 include 'COMMON.GEO'
1440 include 'COMMON.VAR'
1441 include 'COMMON.INTERACT'
1442 include 'COMMON.LOCAL'
1443 include 'COMMON.NAMES'
1444 include 'COMMON.CHAIN'
1445 include 'COMMON.FFIELD'
1446 include 'COMMON.SBRIDGE'
1447 include 'COMMON.HEADER'
1448 include 'COMMON.CONTROL'
1449 include 'COMMON.DBASE'
1450 include 'COMMON.THREAD'
1451 include 'COMMON.TIME1'
1452 C Set up variable list.
1458 if (itype(i).ne.10) then
1460 ialph(i,1)=nvar+nside
1464 if (indphi.gt.0) then
1466 else if (indback.gt.0) then
1471 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1474 c----------------------------------------------------------------------------
1475 subroutine gen_dist_constr
1476 C Generate CA distance constraints.
1477 implicit real*8 (a-h,o-z)
1478 include 'DIMENSIONS'
1479 include 'COMMON.IOUNITS'
1480 include 'COMMON.GEO'
1481 include 'COMMON.VAR'
1482 include 'COMMON.INTERACT'
1483 include 'COMMON.LOCAL'
1484 include 'COMMON.NAMES'
1485 include 'COMMON.CHAIN'
1486 include 'COMMON.FFIELD'
1487 include 'COMMON.SBRIDGE'
1488 include 'COMMON.HEADER'
1489 include 'COMMON.CONTROL'
1490 include 'COMMON.DBASE'
1491 include 'COMMON.THREAD'
1492 include 'COMMON.TIME1'
1493 dimension itype_pdb(maxres)
1494 common /pizda/ itype_pdb
1496 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1497 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1498 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1500 do i=nstart_sup,nstart_sup+nsup-1
1501 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1502 cd & ' seq_pdb', restyp(itype_pdb(i))
1503 do j=i+2,nstart_sup+nsup-1
1505 ihpb(nhpb)=i+nstart_seq-nstart_sup
1506 jhpb(nhpb)=j+nstart_seq-nstart_sup
1508 dhpb(nhpb)=dist(i,j)
1511 cd write (iout,'(a)') 'Distance constraints:'
1516 cd if (ii.gt.nres) then
1521 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1522 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1523 cd & dhpb(i),forcon(i)
1527 c----------------------------------------------------------------------------
1529 implicit real*8 (a-h,o-z)
1530 include 'DIMENSIONS'
1531 include 'COMMON.MAP'
1532 include 'COMMON.IOUNITS'
1533 character*3 angid(4) /'THE','PHI','ALP','OME'/
1534 character*80 mapcard,ucase
1536 read (inp,'(a)') mapcard
1537 mapcard=ucase(mapcard)
1538 if (index(mapcard,'PHI').gt.0) then
1540 else if (index(mapcard,'THE').gt.0) then
1542 else if (index(mapcard,'ALP').gt.0) then
1544 else if (index(mapcard,'OME').gt.0) then
1547 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1548 stop 'Error - illegal variable spec in MAP card.'
1550 call readi (mapcard,'RES1',res1(imap),0)
1551 call readi (mapcard,'RES2',res2(imap),0)
1552 if (res1(imap).eq.0) then
1553 res1(imap)=res2(imap)
1554 else if (res2(imap).eq.0) then
1555 res2(imap)=res1(imap)
1557 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1559 & 'Error - illegal definition of variable group in MAP.'
1560 stop 'Error - illegal definition of variable group in MAP.'
1562 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1563 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1564 call readi(mapcard,'NSTEP',nstep(imap),0)
1565 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1567 & 'Illegal boundary and/or step size specification in MAP.'
1568 stop 'Illegal boundary and/or step size specification in MAP.'
1573 c----------------------------------------------------------------------------
1574 csa subroutine csaread
1575 csa implicit real*8 (a-h,o-z)
1576 csa include 'DIMENSIONS'
1577 csa include 'COMMON.IOUNITS'
1578 csa include 'COMMON.GEO'
1579 csa include 'COMMON.CSA'
1580 csa include 'COMMON.BANK'
1581 csa include 'COMMON.CONTROL'
1582 csa character*80 ucase
1583 csa character*620 mcmcard
1584 csa call card_concat(mcmcard)
1586 csa call readi(mcmcard,'NCONF',nconf,50)
1587 csa call readi(mcmcard,'NADD',nadd,0)
1588 csa call readi(mcmcard,'JSTART',jstart,1)
1589 csa call readi(mcmcard,'JEND',jend,1)
1590 csa call readi(mcmcard,'NSTMAX',nstmax,500000)
1591 csa call readi(mcmcard,'N0',n0,1)
1592 csa call readi(mcmcard,'N1',n1,6)
1593 csa call readi(mcmcard,'N2',n2,4)
1594 csa call readi(mcmcard,'N3',n3,0)
1595 csa call readi(mcmcard,'N4',n4,0)
1596 csa call readi(mcmcard,'N5',n5,0)
1597 csa call readi(mcmcard,'N6',n6,10)
1598 csa call readi(mcmcard,'N7',n7,0)
1599 csa call readi(mcmcard,'N8',n8,0)
1600 csa call readi(mcmcard,'N9',n9,0)
1601 csa call readi(mcmcard,'N14',n14,0)
1602 csa call readi(mcmcard,'N15',n15,0)
1603 csa call readi(mcmcard,'N16',n16,0)
1604 csa call readi(mcmcard,'N17',n17,0)
1605 csa call readi(mcmcard,'N18',n18,0)
1607 csa vdisulf=(index(mcmcard,'DYNSS').gt.0)
1609 csa call readi(mcmcard,'NDIFF',ndiff,2)
1610 csa call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1611 csa call readi(mcmcard,'IS1',is1,1)
1612 csa call readi(mcmcard,'IS2',is2,8)
1613 csa call readi(mcmcard,'NRAN0',nran0,4)
1614 csa call readi(mcmcard,'NRAN1',nran1,2)
1615 csa call readi(mcmcard,'IRR',irr,1)
1616 csa call readi(mcmcard,'NSEED',nseed,20)
1617 csa call readi(mcmcard,'NTOTAL',ntotal,10000)
1618 csa call reada(mcmcard,'CUT1',cut1,2.0d0)
1619 csa call reada(mcmcard,'CUT2',cut2,5.0d0)
1620 csa call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1621 csa call readi(mcmcard,'ICMAX',icmax,3)
1622 csa call readi(mcmcard,'IRESTART',irestart,0)
1623 csac!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1626 csa call reada(mcmcard,'DELE',dele,20.0d0)
1627 csa call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1628 csa call readi(mcmcard,'IREF',iref,0)
1629 csa call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1630 csa call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1631 csa call readi(mcmcard,'NCONF_IN',nconf_in,0)
1632 csa call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1633 csa write (iout,*) "NCONF_IN",nconf_in
1636 c----------------------------------------------------------------------------
1637 cfmc subroutine mcmfread
1638 cfmc implicit real*8 (a-h,o-z)
1639 cfmc include 'DIMENSIONS'
1640 cfmc include 'COMMON.MCMF'
1641 cfmc include 'COMMON.IOUNITS'
1642 cfmc include 'COMMON.GEO'
1643 cfmc character*80 ucase
1644 cfmc character*620 mcmcard
1645 cfmc call card_concat(mcmcard)
1647 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1648 cfmc write(iout,*)'MAXRANT=',maxrant
1649 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1650 cfmc write(iout,*)'MAXFAM=',maxfam
1651 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1652 cfmc write(iout,*)'NNET1=',nnet1
1653 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1654 cfmc write(iout,*)'NNET2=',nnet2
1655 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1656 cfmc write(iout,*)'NNET3=',nnet3
1657 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1658 cfmc write(iout,*)'ILASTT=',ilastt
1659 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1660 cfmc write(iout,*)'MAXSTR=',maxstr
1661 cfmc maxstr_f=maxstr/maxfam
1662 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1663 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1664 cfmc write(iout,*)'NMCMF=',nmcmf
1665 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1666 cfmc write(iout,*)'IFOCUS=',ifocus
1667 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1668 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1669 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1670 cfmc write(iout,*)'INTPRT=',intprt
1671 cfmc call readi(mcmcard,'IPRT',iprt,100)
1672 cfmc write(iout,*)'IPRT=',iprt
1673 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1674 cfmc write(iout,*)'IMAXTR=',imaxtr
1675 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1676 cfmc write(iout,*)'MAXEVEN=',maxeven
1677 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1678 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1679 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1680 cfmc write(iout,*)'INIMIN=',inimin
1681 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1682 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1683 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1684 cfmc write(iout,*)'NTHREAD=',nthread
1685 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1686 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1687 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1688 cfmc write(iout,*)'MAXPERT=',maxpert
1689 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1690 cfmc write(iout,*)'IRMSD=',irmsd
1691 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1692 cfmc write(iout,*)'DENEMIN=',denemin
1693 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1694 cfmc write(iout,*)'RCUT1S=',rcut1s
1695 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1696 cfmc write(iout,*)'RCUT1E=',rcut1e
1697 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1698 cfmc write(iout,*)'RCUT2S=',rcut2s
1699 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1700 cfmc write(iout,*)'RCUT2E=',rcut2e
1701 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1702 cfmc write(iout,*)'DPERT1=',d_pert1
1703 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1704 cfmc write(iout,*)'DPERT1A=',d_pert1a
1705 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1706 cfmc write(iout,*)'DPERT2=',d_pert2
1707 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1708 cfmc write(iout,*)'DPERT2A=',d_pert2a
1709 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1710 cfmc write(iout,*)'DPERT2B=',d_pert2b
1711 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1712 cfmc write(iout,*)'DPERT2C=',d_pert2c
1713 cfmc d_pert1=deg2rad*d_pert1
1714 cfmc d_pert1a=deg2rad*d_pert1a
1715 cfmc d_pert2=deg2rad*d_pert2
1716 cfmc d_pert2a=deg2rad*d_pert2a
1717 cfmc d_pert2b=deg2rad*d_pert2b
1718 cfmc d_pert2c=deg2rad*d_pert2c
1719 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1720 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1721 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1722 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1723 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1724 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1725 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1726 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1727 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1728 cfmc write(iout,*)'RCUTINI=',rcutini
1729 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1730 cfmc write(iout,*)'GRAT=',grat
1731 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1732 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1736 c----------------------------------------------------------------------------
1738 implicit real*8 (a-h,o-z)
1739 include 'DIMENSIONS'
1740 include 'COMMON.MCM'
1741 include 'COMMON.MCE'
1742 include 'COMMON.IOUNITS'
1744 character*320 mcmcard
1745 call card_concat(mcmcard)
1746 call readi(mcmcard,'MAXACC',maxacc,100)
1747 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1748 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1749 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1750 call readi(mcmcard,'MAXREPM',maxrepm,200)
1751 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1752 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1753 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1754 call reada(mcmcard,'E_UP',e_up,5.0D0)
1755 call reada(mcmcard,'DELTE',delte,0.1D0)
1756 call readi(mcmcard,'NSWEEP',nsweep,5)
1757 call readi(mcmcard,'NSTEPH',nsteph,0)
1758 call readi(mcmcard,'NSTEPC',nstepc,0)
1759 call reada(mcmcard,'TMIN',tmin,298.0D0)
1760 call reada(mcmcard,'TMAX',tmax,298.0D0)
1761 call readi(mcmcard,'NWINDOW',nwindow,0)
1762 call readi(mcmcard,'PRINT_MC',print_mc,0)
1763 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1764 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1765 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1766 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1767 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1768 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1769 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1770 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1771 if (nwindow.gt.0) then
1772 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1774 winlen(i)=winend(i)-winstart(i)+1
1777 if (tmax.lt.tmin) tmax=tmin
1778 if (tmax.eq.tmin) then
1782 if (nstepc.gt.0 .and. nsteph.gt.0) then
1783 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1784 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1786 C Probabilities of different move types
1787 sumpro_type(0)=0.0D0
1788 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1789 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1790 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1791 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1792 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1793 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1794 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1796 print *,'i',i,' sumprotype',sumpro_type(i)
1797 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1798 print *,'i',i,' sumprotype',sumpro_type(i)
1802 c----------------------------------------------------------------------------
1803 subroutine read_minim
1804 implicit real*8 (a-h,o-z)
1805 include 'DIMENSIONS'
1806 include 'COMMON.MINIM'
1807 include 'COMMON.IOUNITS'
1809 character*320 minimcard
1810 call card_concat(minimcard)
1811 call readi(minimcard,'MAXMIN',maxmin,2000)
1812 call readi(minimcard,'MAXFUN',maxfun,5000)
1813 call readi(minimcard,'MINMIN',minmin,maxmin)
1814 call readi(minimcard,'MINFUN',minfun,maxmin)
1815 call reada(minimcard,'TOLF',tolf,1.0D-2)
1816 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1817 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1818 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1819 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1820 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1821 & 'Options in energy minimization:'
1822 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1823 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1824 & 'MinMin:',MinMin,' MinFun:',MinFun,
1825 & ' TolF:',TolF,' RTolF:',RTolF
1828 c----------------------------------------------------------------------------
1829 subroutine read_angles(kanal,*)
1830 implicit real*8 (a-h,o-z)
1831 include 'DIMENSIONS'
1832 include 'COMMON.GEO'
1833 include 'COMMON.VAR'
1834 include 'COMMON.CHAIN'
1835 include 'COMMON.IOUNITS'
1836 include 'COMMON.CONTROL'
1837 c Read angles from input
1839 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1840 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1841 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1842 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1845 c 9/7/01 avoid 180 deg valence angle
1846 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1848 theta(i)=deg2rad*theta(i)
1849 phi(i)=deg2rad*phi(i)
1850 alph(i)=deg2rad*alph(i)
1851 omeg(i)=deg2rad*omeg(i)
1856 c----------------------------------------------------------------------------
1857 subroutine reada(rekord,lancuch,wartosc,default)
1859 character*(*) rekord,lancuch
1860 double precision wartosc,default
1863 iread=index(rekord,lancuch)
1864 if (iread.eq.0) then
1868 iread=iread+ilen(lancuch)+1
1869 read (rekord(iread:),*,err=10,end=10) wartosc
1874 c----------------------------------------------------------------------------
1875 subroutine readi(rekord,lancuch,wartosc,default)
1877 character*(*) rekord,lancuch
1878 integer wartosc,default
1881 iread=index(rekord,lancuch)
1882 if (iread.eq.0) then
1886 iread=iread+ilen(lancuch)+1
1887 read (rekord(iread:),*,err=10,end=10) wartosc
1892 c----------------------------------------------------------------------------
1893 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1896 integer tablica(dim),default
1897 character*(*) rekord,lancuch
1904 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1905 if (iread.eq.0) return
1906 iread=iread+ilen(lancuch)+1
1907 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1910 c----------------------------------------------------------------------------
1911 subroutine multreada(rekord,lancuch,tablica,dim,default)
1914 double precision tablica(dim),default
1915 character*(*) rekord,lancuch
1922 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1923 if (iread.eq.0) return
1924 iread=iread+ilen(lancuch)+1
1925 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1928 c----------------------------------------------------------------------------
1929 subroutine openunits
1930 implicit real*8 (a-h,o-z)
1931 include 'DIMENSIONS'
1934 character*16 form,nodename
1937 include 'COMMON.SETUP'
1938 include 'COMMON.IOUNITS'
1940 include 'COMMON.CONTROL'
1941 integer lenpre,lenpot,ilen,lentmp
1943 character*3 out1file_text,ucase
1946 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1947 call getenv_loc("PREFIX",prefix)
1949 call getenv_loc("POT",pot)
1950 call getenv_loc("DIRTMP",tmpdir)
1951 call getenv_loc("CURDIR",curdir)
1952 call getenv_loc("OUT1FILE",out1file_text)
1953 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1954 out1file_text=ucase(out1file_text)
1955 if (out1file_text(1:1).eq."Y") then
1958 out1file=fg_rank.gt.0
1963 if (lentmp.gt.0) then
1964 write (*,'(80(1h!))')
1965 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1966 write (*,'(80(1h!))')
1967 write (*,*)"All output files will be on node /tmp directory."
1969 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1970 if (me.eq.king) then
1971 write (*,*) "The master node is ",nodename
1972 else if (fg_rank.eq.0) then
1973 write (*,*) "I am the CG slave node ",nodename
1975 write (*,*) "I am the FG slave node ",nodename
1978 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1979 lenpre = lentmp+lenpre+1
1981 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1982 C Get the names and open the input files
1983 #if defined(WINIFL) || defined(WINPGI)
1984 open(1,file=pref_orig(:ilen(pref_orig))//
1985 & '.inp',status='old',readonly,shared)
1986 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1987 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1988 C Get parameter filenames and open the parameter files.
1989 call getenv_loc('BONDPAR',bondname)
1990 open (ibond,file=bondname,status='old',readonly,shared)
1991 call getenv_loc('THETPAR',thetname)
1992 open (ithep,file=thetname,status='old',readonly,shared)
1994 call getenv_loc('THETPARPDB',thetname_pdb)
1995 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1997 call getenv_loc('ROTPAR',rotname)
1998 open (irotam,file=rotname,status='old',readonly,shared)
2000 call getenv_loc('ROTPARPDB',rotname_pdb)
2001 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2003 call getenv_loc('TORPAR',torname)
2004 open (itorp,file=torname,status='old',readonly,shared)
2005 call getenv_loc('TORDPAR',tordname)
2006 open (itordp,file=tordname,status='old',readonly,shared)
2007 call getenv_loc('FOURIER',fouriername)
2008 open (ifourier,file=fouriername,status='old',readonly,shared)
2009 call getenv_loc('ELEPAR',elename)
2010 open (ielep,file=elename,status='old',readonly,shared)
2011 call getenv_loc('SIDEPAR',sidename)
2012 open (isidep,file=sidename,status='old',readonly,shared)
2013 #elif (defined CRAY) || (defined AIX)
2014 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2016 c print *,"Processor",myrank," opened file 1"
2017 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2018 c print *,"Processor",myrank," opened file 9"
2019 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2020 C Get parameter filenames and open the parameter files.
2021 call getenv_loc('BONDPAR',bondname)
2022 open (ibond,file=bondname,status='old',action='read')
2023 c print *,"Processor",myrank," opened file IBOND"
2024 call getenv_loc('THETPAR',thetname)
2025 open (ithep,file=thetname,status='old',action='read')
2026 c print *,"Processor",myrank," opened file ITHEP"
2028 call getenv_loc('THETPARPDB',thetname_pdb)
2029 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2031 call getenv_loc('ROTPAR',rotname)
2032 open (irotam,file=rotname,status='old',action='read')
2033 c print *,"Processor",myrank," opened file IROTAM"
2035 call getenv_loc('ROTPARPDB',rotname_pdb)
2036 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2038 call getenv_loc('TORPAR',torname)
2039 open (itorp,file=torname,status='old',action='read')
2040 c print *,"Processor",myrank," opened file ITORP"
2041 call getenv_loc('TORDPAR',tordname)
2042 open (itordp,file=tordname,status='old',action='read')
2043 c print *,"Processor",myrank," opened file ITORDP"
2044 call getenv_loc('SCCORPAR',sccorname)
2045 open (isccor,file=sccorname,status='old',action='read')
2046 c print *,"Processor",myrank," opened file ISCCOR"
2047 call getenv_loc('FOURIER',fouriername)
2048 open (ifourier,file=fouriername,status='old',action='read')
2049 c print *,"Processor",myrank," opened file IFOURIER"
2050 call getenv_loc('ELEPAR',elename)
2051 open (ielep,file=elename,status='old',action='read')
2052 c print *,"Processor",myrank," opened file IELEP"
2053 call getenv_loc('SIDEPAR',sidename)
2054 open (isidep,file=sidename,status='old',action='read')
2055 c print *,"Processor",myrank," opened file ISIDEP"
2056 c print *,"Processor",myrank," opened parameter files"
2058 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2059 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2060 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2061 C Get parameter filenames and open the parameter files.
2062 call getenv_loc('BONDPAR',bondname)
2063 open (ibond,file=bondname,status='old')
2064 call getenv_loc('THETPAR',thetname)
2065 open (ithep,file=thetname,status='old')
2067 call getenv_loc('THETPARPDB',thetname_pdb)
2068 open (ithep_pdb,file=thetname_pdb,status='old')
2070 call getenv_loc('ROTPAR',rotname)
2071 open (irotam,file=rotname,status='old')
2073 call getenv_loc('ROTPARPDB',rotname_pdb)
2074 open (irotam_pdb,file=rotname_pdb,status='old')
2076 call getenv_loc('TORPAR',torname)
2077 open (itorp,file=torname,status='old')
2078 call getenv_loc('TORDPAR',tordname)
2079 open (itordp,file=tordname,status='old')
2080 call getenv_loc('SCCORPAR',sccorname)
2081 open (isccor,file=sccorname,status='old')
2082 call getenv_loc('FOURIER',fouriername)
2083 open (ifourier,file=fouriername,status='old')
2084 call getenv_loc('ELEPAR',elename)
2085 open (ielep,file=elename,status='old')
2086 call getenv_loc('SIDEPAR',sidename)
2087 open (isidep,file=sidename,status='old')
2089 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2091 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2092 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2093 C Get parameter filenames and open the parameter files.
2094 call getenv_loc('BONDPAR',bondname)
2095 open (ibond,file=bondname,status='old',action='read')
2096 call getenv_loc('THETPAR',thetname)
2097 open (ithep,file=thetname,status='old',action='read')
2099 call getenv_loc('THETPARPDB',thetname_pdb)
2100 print *,"thetname_pdb ",thetname_pdb
2101 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2102 print *,ithep_pdb," opened"
2104 call getenv_loc('ROTPAR',rotname)
2105 open (irotam,file=rotname,status='old',action='read')
2107 call getenv_loc('ROTPARPDB',rotname_pdb)
2108 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2110 call getenv_loc('TORPAR',torname)
2111 open (itorp,file=torname,status='old',action='read')
2112 call getenv_loc('TORDPAR',tordname)
2113 open (itordp,file=tordname,status='old',action='read')
2114 call getenv_loc('SCCORPAR',sccorname)
2115 open (isccor,file=sccorname,status='old',action='read')
2116 call getenv_loc('FOURIER',fouriername)
2117 open (ifourier,file=fouriername,status='old',action='read')
2118 call getenv_loc('ELEPAR',elename)
2119 open (ielep,file=elename,status='old',action='read')
2120 call getenv_loc('SIDEPAR',sidename)
2121 open (isidep,file=sidename,status='old',action='read')
2125 C 8/9/01 In the newest version SCp interaction constants are read from a file
2126 C Use -DOLDSCP to use hard-coded constants instead.
2128 call getenv_loc('SCPPAR',scpname)
2129 #if defined(WINIFL) || defined(WINPGI)
2130 open (iscpp,file=scpname,status='old',readonly,shared)
2131 #elif (defined CRAY) || (defined AIX)
2132 open (iscpp,file=scpname,status='old',action='read')
2134 open (iscpp,file=scpname,status='old')
2136 open (iscpp,file=scpname,status='old',action='read')
2139 call getenv_loc('PATTERN',patname)
2140 #if defined(WINIFL) || defined(WINPGI)
2141 open (icbase,file=patname,status='old',readonly,shared)
2142 #elif (defined CRAY) || (defined AIX)
2143 open (icbase,file=patname,status='old',action='read')
2145 open (icbase,file=patname,status='old')
2147 open (icbase,file=patname,status='old',action='read')
2150 C Open output file only for CG processes
2151 c print *,"Processor",myrank," fg_rank",fg_rank
2152 if (fg_rank.eq.0) then
2154 if (nodes.eq.1) then
2157 npos = dlog10(dfloat(nodes-1))+1
2159 if (npos.lt.3) npos=3
2160 write (liczba,'(i1)') npos
2161 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2163 write (liczba,form) me
2164 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2165 & liczba(:ilen(liczba))
2166 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2168 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2170 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2171 & liczba(:ilen(liczba))//'.mol2'
2172 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2173 & liczba(:ilen(liczba))//'.stat'
2175 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2176 & //liczba(:ilen(liczba))//'.stat')
2177 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2180 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2181 & liczba(:ilen(liczba))//'.const'
2186 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2187 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2188 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2189 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2190 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2192 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2194 rest2name=prefix(:ilen(prefix))//'.rst'
2196 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2199 #if defined(AIX) || defined(PGI)
2200 if (me.eq.king .or. .not. out1file)
2201 & open(iout,file=outname,status='unknown')
2204 if (fg_rank.gt.0) then
2205 write (liczba,'(i3.3)') myrank/nfgtasks
2206 write (ll,'(bz,i3.3)') fg_rank
2207 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2213 open(igeom,file=intname,status='unknown',position='append')
2214 open(ipdb,file=pdbname,status='unknown')
2215 open(imol2,file=mol2name,status='unknown')
2216 open(istat,file=statname,status='unknown',position='append')
2218 c1out open(iout,file=outname,status='unknown')
2221 if (me.eq.king .or. .not.out1file)
2222 & open(iout,file=outname,status='unknown')
2225 if (fg_rank.gt.0) then
2226 print "Processor",fg_rank," opening output file"
2227 write (liczba,'(i3.3)') myrank/nfgtasks
2228 write (ll,'(bz,i3.3)') fg_rank
2229 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2235 open(igeom,file=intname,status='unknown',access='append')
2236 open(ipdb,file=pdbname,status='unknown')
2237 open(imol2,file=mol2name,status='unknown')
2238 open(istat,file=statname,status='unknown',access='append')
2240 c1out open(iout,file=outname,status='unknown')
2243 csa csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2244 csa csa_seed=prefix(:lenpre)//'.CSA.seed'
2245 csa csa_history=prefix(:lenpre)//'.CSA.history'
2246 csa csa_bank=prefix(:lenpre)//'.CSA.bank'
2247 csa csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2248 csa csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2249 csa csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2250 csac!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2251 csa csa_int=prefix(:lenpre)//'.int'
2252 csa csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2253 csa csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2254 csa csa_in=prefix(:lenpre)//'.CSA.in'
2255 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2258 write (iout,'(80(1h-))')
2259 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2260 write (iout,'(80(1h-))')
2261 write (iout,*) "Input file : ",
2262 & pref_orig(:ilen(pref_orig))//'.inp'
2263 write (iout,*) "Output file : ",
2264 & outname(:ilen(outname))
2266 write (iout,*) "Sidechain potential file : ",
2267 & sidename(:ilen(sidename))
2269 write (iout,*) "SCp potential file : ",
2270 & scpname(:ilen(scpname))
2272 write (iout,*) "Electrostatic potential file : ",
2273 & elename(:ilen(elename))
2274 write (iout,*) "Cumulant coefficient file : ",
2275 & fouriername(:ilen(fouriername))
2276 write (iout,*) "Torsional parameter file : ",
2277 & torname(:ilen(torname))
2278 write (iout,*) "Double torsional parameter file : ",
2279 & tordname(:ilen(tordname))
2280 write (iout,*) "SCCOR parameter file : ",
2281 & sccorname(:ilen(sccorname))
2282 write (iout,*) "Bond & inertia constant file : ",
2283 & bondname(:ilen(bondname))
2284 write (iout,*) "Bending parameter file : ",
2285 & thetname(:ilen(thetname))
2286 write (iout,*) "Rotamer parameter file : ",
2287 & rotname(:ilen(rotname))
2288 write (iout,*) "Threading database : ",
2289 & patname(:ilen(patname))
2291 &write (iout,*)" DIRTMP : ",
2293 write (iout,'(80(1h-))')
2297 c----------------------------------------------------------------------------
2298 subroutine card_concat(card)
2299 implicit real*8 (a-h,o-z)
2300 include 'DIMENSIONS'
2301 include 'COMMON.IOUNITS'
2303 character*80 karta,ucase
2305 read (inp,'(a)') karta
2308 do while (karta(80:80).eq.'&')
2309 card=card(:ilen(card)+1)//karta(:79)
2310 read (inp,'(a)') karta
2313 card=card(:ilen(card)+1)//karta
2316 c----------------------------------------------------------------------------------
2318 implicit real*8 (a-h,o-z)
2319 include 'DIMENSIONS'
2320 include 'COMMON.CHAIN'
2321 include 'COMMON.IOUNITS'
2323 open(irest2,file=rest2name,status='unknown')
2324 read(irest2,*) totT,EK,potE,totE,t_bath
2326 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2329 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2332 read (irest2,*) iset
2337 c---------------------------------------------------------------------------------
2338 subroutine read_fragments
2339 implicit real*8 (a-h,o-z)
2340 include 'DIMENSIONS'
2344 include 'COMMON.SETUP'
2345 include 'COMMON.CHAIN'
2346 include 'COMMON.IOUNITS'
2348 include 'COMMON.CONTROL'
2349 read(inp,*) nset,nfrag,npair,nfrag_back
2350 if(me.eq.king.or..not.out1file)
2351 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2352 & " nfrag_back",nfrag_back
2354 read(inp,*) mset(iset)
2356 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2358 if(me.eq.king.or..not.out1file)
2359 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2360 & ifrag(2,i,iset), qinfrag(i,iset)
2363 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2365 if(me.eq.king.or..not.out1file)
2366 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2367 & ipair(2,i,iset), qinpair(i,iset)
2370 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2371 & wfrag_back(3,i,iset),
2372 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2373 if(me.eq.king.or..not.out1file)
2374 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2375 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2380 c-------------------------------------------------------------------------------
2381 subroutine read_dist_constr
2382 implicit real*8 (a-h,o-z)
2383 include 'DIMENSIONS'
2387 include 'COMMON.SETUP'
2388 include 'COMMON.CONTROL'
2389 include 'COMMON.CHAIN'
2390 include 'COMMON.IOUNITS'
2391 include 'COMMON.SBRIDGE'
2392 integer ifrag_(2,100),ipair_(2,100)
2393 double precision wfrag_(100),wpair_(100)
2394 character*500 controlcard
2395 c write (iout,*) "Calling read_dist_constr"
2396 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2398 call card_concat(controlcard)
2399 call readi(controlcard,"NFRAG",nfrag_,0)
2400 call readi(controlcard,"NPAIR",npair_,0)
2401 call readi(controlcard,"NDIST",ndist_,0)
2402 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2403 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2404 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2405 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2406 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2407 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2408 c write (iout,*) "IFRAG"
2410 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2412 c write (iout,*) "IPAIR"
2414 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2416 if (.not.refstr .and. nfrag.gt.0) then
2418 & "ERROR: no reference structure to compute distance restraints"
2420 & "Restraints must be specified explicitly (NDIST=number)"
2423 if (nfrag.lt.2 .and. npair.gt.0) then
2424 write (iout,*) "ERROR: Less than 2 fragments specified",
2425 & " but distance restraints between pairs requested"
2430 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2431 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2432 & ifrag_(2,i)=nstart_sup+nsup-1
2433 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2435 if (wfrag_(i).gt.0.0d0) then
2436 do j=ifrag_(1,i),ifrag_(2,i)-1
2437 do k=j+1,ifrag_(2,i)
2438 write (iout,*) "j",j," k",k
2440 if (constr_dist.eq.1) then
2445 forcon(nhpb)=wfrag_(i)
2446 else if (constr_dist.eq.2) then
2447 if (ddjk.le.dist_cut) then
2452 forcon(nhpb)=wfrag_(i)
2459 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2462 if (.not.out1file .or. me.eq.king)
2463 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2464 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2466 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2467 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2474 if (wpair_(i).gt.0.0d0) then
2482 do j=ifrag_(1,ii),ifrag_(2,ii)
2483 do k=ifrag_(1,jj),ifrag_(2,jj)
2487 forcon(nhpb)=wpair_(i)
2488 dhpb(nhpb)=dist(j,k)
2490 if (.not.out1file .or. me.eq.king)
2491 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2492 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2494 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2495 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2502 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2503 & ibecarb(i),forcon(nhpb+1)
2504 if (forcon(nhpb+1).gt.0.0d0) then
2506 if (ibecarb(i).gt.0) then
2507 ihpb(i)=ihpb(i)+nres
2508 jhpb(i)=jhpb(i)+nres
2510 if (dhpb(nhpb).eq.0.0d0)
2511 & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2515 if (.not.out1file .or. me.eq.king) then
2518 write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2519 & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2527 c-------------------------------------------------------------------------------
2529 subroutine flush(iu)
2534 subroutine flush(iu)
2539 c------------------------------------------------------------------------------
2540 subroutine copy_to_tmp(source)
2541 include "DIMENSIONS"
2542 include "COMMON.IOUNITS"
2543 character*(*) source
2544 character* 256 tmpfile
2548 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2549 inquire(file=tmpfile,exist=ex)
2551 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2552 & " to temporary directory..."
2553 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2554 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2558 c------------------------------------------------------------------------------
2559 subroutine move_from_tmp(source)
2560 include "DIMENSIONS"
2561 include "COMMON.IOUNITS"
2562 character*(*) source
2565 write (*,*) "Moving ",source(:ilen(source)),
2566 & " from temporary directory to working directory"
2567 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2568 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2571 c------------------------------------------------------------------------------
2572 subroutine random_init(seed)
2574 C Initialize random number generator
2576 implicit real*8 (a-h,o-z)
2577 include 'DIMENSIONS'
2583 logical OKRandom, prng_restart
2585 integer iseed_array(4)
2587 include 'COMMON.IOUNITS'
2588 include 'COMMON.TIME1'
2589 include 'COMMON.THREAD'
2590 include 'COMMON.SBRIDGE'
2591 include 'COMMON.CONTROL'
2592 include 'COMMON.MCM'
2593 include 'COMMON.MAP'
2594 include 'COMMON.HEADER'
2595 csa include 'COMMON.CSA'
2596 include 'COMMON.CHAIN'
2597 include 'COMMON.MUCA'
2599 include 'COMMON.FFIELD'
2600 include 'COMMON.SETUP'
2601 iseed=-dint(dabs(seed))
2602 if (iseed.eq.0) then
2603 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2604 & 'Random seed undefined. The program will stop.'
2605 write (*,'(/80(1h*)/20x,a/80(1h*))')
2606 & 'Random seed undefined. The program will stop.'
2608 call mpi_finalize(mpi_comm_world,ierr)
2610 stop 'Bad random seed.'
2613 if (fg_rank.eq.0) then
2617 if(me.eq.king .or. .not. out1file)
2618 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2619 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2620 OKRandom = prng_restart(me,iseedi8)
2623 tmp=65536.0d0**(4-i)
2624 iseed_array(i) = dint(seed/tmp)
2625 seed=seed-iseed_array(i)*tmp
2627 if(me.eq.king .or. .not. out1file)
2628 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2629 & (iseed_array(i),i=1,4)
2630 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2631 & (iseed_array(i),i=1,4)
2632 OKRandom = prng_restart(me,iseed_array)
2635 r1=ran_number(0.0D0,1.0D0)
2636 if(me.eq.king .or. .not. out1file)
2637 & write (iout,*) 'ran_num',r1
2638 if (r1.lt.0.0d0) OKRandom=.false.
2640 if (.not.OKRandom) then
2641 write (iout,*) 'PRNG IS NOT WORKING!!!'
2642 print *,'PRNG IS NOT WORKING!!!'
2645 call mpi_abort(mpi_comm_world,error_msg,ierr)
2648 write (iout,*) 'too many processors for parallel prng'
2649 write (*,*) 'too many processors for parallel prng'
2657 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)