2 implicit real*8 (a-h,o-z)
8 include 'COMMON.CONTROL'
9 include 'COMMON.SBRIDGE'
10 include 'COMMON.IOUNITS'
11 include 'COMMON.SPLITELE'
13 C Read force-field parameters except weights
15 C Read job setup parameters
17 C Read control parameters for energy minimzation if required
18 if (minim) call read_minim
19 C Read MCM control parameters if required
20 if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
21 C Read MD control parameters if reqjuired
22 if (modecalc.eq.12) call read_MDpar
23 C Read MREMD control parameters if required
24 if (modecalc.eq.14) then
28 C Read MUCA control parameters if required
29 if (lmuca) call read_muca
30 C Read CSA control parameters if required (from fort.40 if exists
31 C otherwise from general input file)
32 if (modecalc.eq.8) then
33 inquire (file="fort.40",exist=file_exist)
34 if (.not.file_exist) call csaread
36 cfmc if (modecalc.eq.10) call mcmfread
37 C Read molecule information, molecule geometry, energy-term weights, and
38 C restraints if requested
40 C Print restraint information
42 if (.not. out1file .or. me.eq.king) then
45 &write (iout,'(a,i5,a)') "The following",nhpb-nss,
46 & " distance constraints have been imposed"
48 write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(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'
77 include 'COMMON.CHAIN'
80 include 'COMMON.FFIELD'
81 include 'COMMON.INTERACT'
82 include 'COMMON.SETUP'
83 include 'COMMON.SPLITELE'
84 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
85 character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
87 character*320 controlcard
92 read (INP,'(a)') titel
93 call card_concat(controlcard)
94 c out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
95 c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
96 call reada(controlcard,'SEED',seed,0.0D0)
97 call random_init(seed)
98 C Set up the time limit (caution! The time must be input in minutes!)
99 read_cart=index(controlcard,'READ_CART').gt.0
100 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
101 call readi(controlcard,'SYM',symetr,1)
102 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
103 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
104 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
105 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
106 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
107 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
108 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
109 call reada(controlcard,'DRMS',drms,0.1D0)
110 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
111 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
112 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
113 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
114 write (iout,'(a,f10.1)')'DRMS = ',drms
115 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
116 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
118 call readi(controlcard,'NZ_START',nz_start,0)
119 call readi(controlcard,'NZ_END',nz_end,0)
120 call readi(controlcard,'IZ_SC',iz_sc,0)
122 safety = 60.0d0*safety
125 call reada(controlcard,"T_BATH",t_bath,300.0d0)
126 minim=(index(controlcard,'MINIMIZE').gt.0)
127 dccart=(index(controlcard,'CART').gt.0)
128 overlapsc=(index(controlcard,'OVERLAP').gt.0)
129 overlapsc=.not.overlapsc
130 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
131 searchsc=.not.searchsc
132 sideadd=(index(controlcard,'SIDEADD').gt.0)
133 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
134 outpdb=(index(controlcard,'PDBOUT').gt.0)
135 outmol2=(index(controlcard,'MOL2OUT').gt.0)
136 pdbref=(index(controlcard,'PDBREF').gt.0)
137 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
138 indpdb=index(controlcard,'PDBSTART')
139 extconf=(index(controlcard,'EXTCONF').gt.0)
140 call readi(controlcard,'IPRINT',iprint,0)
141 call readi(controlcard,'MAXGEN',maxgen,10000)
142 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
143 call readi(controlcard,"KDIAG",kdiag,0)
144 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
145 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
146 & write (iout,*) "RESCALE_MODE",rescale_mode
147 split_ene=index(controlcard,'SPLIT_ENE').gt.0
148 if (index(controlcard,'REGULAR').gt.0.0D0) then
149 call reada(controlcard,'WEIDIS',weidis,0.1D0)
153 if (index(controlcard,'CHECKGRAD').gt.0) then
155 if (index(controlcard,'CART').gt.0) then
157 elseif (index(controlcard,'CARINT').gt.0) then
162 elseif (index(controlcard,'THREAD').gt.0) then
164 call readi(controlcard,'THREAD',nthread,0)
165 if (nthread.gt.0) then
166 call reada(controlcard,'WEIDIS',weidis,0.1D0)
169 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
170 stop 'Error termination in Read_Control.'
172 else if (index(controlcard,'MCMA').gt.0) then
174 else if (index(controlcard,'MCEE').gt.0) then
176 else if (index(controlcard,'MULTCONF').gt.0) then
178 else if (index(controlcard,'MAP').gt.0) then
180 call readi(controlcard,'MAP',nmap,0)
181 else if (index(controlcard,'CSA').gt.0) then
183 crc else if (index(controlcard,'ZSCORE').gt.0) then
185 crc ZSCORE is rm from UNRES, modecalc=9 is available
188 cfcm else if (index(controlcard,'MCMF').gt.0) then
190 else if (index(controlcard,'SOFTREG').gt.0) then
192 else if (index(controlcard,'CHECK_BOND').gt.0) then
194 else if (index(controlcard,'TEST').gt.0) then
196 else if (index(controlcard,'MD').gt.0) then
198 else if (index(controlcard,'RE ').gt.0) then
202 lmuca=index(controlcard,'MUCA').gt.0
203 call readi(controlcard,'MUCADYN',mucadyn,0)
204 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
205 if (lmuca .and. (me.eq.king .or. .not.out1file ))
207 write (iout,*) 'MUCADYN=',mucadyn
208 write (iout,*) 'MUCASMOOTH=',muca_smooth
211 iscode=index(controlcard,'ONE_LETTER')
212 indphi=index(controlcard,'PHI')
213 indback=index(controlcard,'BACK')
214 iranconf=index(controlcard,'RAND_CONF')
215 i2ndstr=index(controlcard,'USE_SEC_PRED')
216 gradout=index(controlcard,'GRADOUT').gt.0
217 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
218 C DISTCHAINMAX become obsolete for periodic boundry condition
219 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
220 C Reading the dimensions of box in x,y,z coordinates
221 call reada(controlcard,'BOXX',boxxsize,100.0d0)
222 call reada(controlcard,'BOXY',boxysize,100.0d0)
223 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
224 c Cutoff range for interactions
225 call reada(controlcard,"R_CUT",r_cut,15.0d0)
226 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
227 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
228 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
229 if (lipthick.gt.0.0d0) then
230 bordliptop=(boxzsize+lipthick)/2.0
231 bordlipbot=bordliptop-lipthick
233 if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0))
234 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
235 buflipbot=bordlipbot+lipbufthick
236 bufliptop=bordliptop-lipbufthick
237 if ((lipbufthick*2.0d0).gt.lipthick)
238 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
240 if (me.eq.king .or. .not.out1file )
241 & write (iout,*) "DISTCHAINMAX",distchainmax
243 if(me.eq.king.or..not.out1file)
244 & write (iout,'(2a)') diagmeth(kdiag),
245 & ' routine used to diagonalize matrices.'
248 c--------------------------------------------------------------------------
249 subroutine read_REMDpar
253 implicit real*8 (a-h,o-z)
255 include 'COMMON.IOUNITS'
256 include 'COMMON.TIME1'
259 include 'COMMON.LANGEVIN'
261 include 'COMMON.LANGEVIN.lang0'
263 include 'COMMON.INTERACT'
264 include 'COMMON.NAMES'
266 include 'COMMON.REMD'
267 include 'COMMON.CONTROL'
268 include 'COMMON.SETUP'
270 character*320 controlcard
271 character*3200 controlcard1
272 integer iremd_m_total
274 if(me.eq.king.or..not.out1file)
275 & write (iout,*) "REMD setup"
277 call card_concat(controlcard)
278 call readi(controlcard,"NREP",nrep,3)
279 call readi(controlcard,"NSTEX",nstex,1000)
280 call reada(controlcard,"RETMIN",retmin,10.0d0)
281 call reada(controlcard,"RETMAX",retmax,1000.0d0)
282 mremdsync=(index(controlcard,'SYNC').gt.0)
283 call readi(controlcard,"NSYN",i_sync_step,100)
284 restart1file=(index(controlcard,'REST1FILE').gt.0)
285 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
286 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
287 if(max_cache_traj_use.gt.max_cache_traj)
288 & max_cache_traj_use=max_cache_traj
289 if(me.eq.king.or..not.out1file) then
290 cd if (traj1file) then
291 crc caching is in testing - NTWX is not ignored
292 cd write (iout,*) "NTWX value is ignored"
293 cd write (iout,*) " trajectory is stored to one file by master"
294 cd write (iout,*) " before exchange at NSTEX intervals"
296 write (iout,*) "NREP= ",nrep
297 write (iout,*) "NSTEX= ",nstex
298 write (iout,*) "SYNC= ",mremdsync
299 write (iout,*) "NSYN= ",i_sync_step
300 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
303 if (index(controlcard,'TLIST').gt.0) then
305 call card_concat(controlcard1)
306 read(controlcard1,*) (remd_t(i),i=1,nrep)
307 if(me.eq.king.or..not.out1file)
308 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
311 if (index(controlcard,'MLIST').gt.0) then
313 call card_concat(controlcard1)
314 read(controlcard1,*) (remd_m(i),i=1,nrep)
315 if(me.eq.king.or..not.out1file) then
316 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
319 iremd_m_total=iremd_m_total+remd_m(i)
321 write (iout,*) 'Total number of replicas ',iremd_m_total
324 if(me.eq.king.or..not.out1file)
325 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
328 c--------------------------------------------------------------------------
329 subroutine read_MDpar
333 implicit real*8 (a-h,o-z)
335 include 'COMMON.IOUNITS'
336 include 'COMMON.TIME1'
339 include 'COMMON.LANGEVIN'
341 include 'COMMON.LANGEVIN.lang0'
343 include 'COMMON.INTERACT'
344 include 'COMMON.NAMES'
346 include 'COMMON.SETUP'
347 include 'COMMON.CONTROL'
348 include 'COMMON.SPLITELE'
350 character*320 controlcard
352 call card_concat(controlcard)
353 call readi(controlcard,"NSTEP",n_timestep,1000000)
354 call readi(controlcard,"NTWE",ntwe,100)
355 call readi(controlcard,"NTWX",ntwx,1000)
356 call reada(controlcard,"DT",d_time,1.0d-1)
357 call reada(controlcard,"DVMAX",dvmax,2.0d1)
358 call reada(controlcard,"DAMAX",damax,1.0d1)
359 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
360 call readi(controlcard,"LANG",lang,0)
361 RESPA = index(controlcard,"RESPA") .gt. 0
362 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
363 ntime_split0=ntime_split
364 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
365 ntime_split0=ntime_split
366 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
367 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
368 rest = index(controlcard,"REST").gt.0
369 tbf = index(controlcard,"TBF").gt.0
370 usampl = index(controlcard,"USAMPL").gt.0
372 mdpdb = index(controlcard,"MDPDB").gt.0
373 call reada(controlcard,"T_BATH",t_bath,300.0d0)
374 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
375 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
376 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
377 if (count_reset_moment.eq.0) count_reset_moment=1000000000
378 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
379 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
380 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
381 if (count_reset_vel.eq.0) count_reset_vel=1000000000
382 large = index(controlcard,"LARGE").gt.0
383 print_compon = index(controlcard,"PRINT_COMPON").gt.0
384 rattle = index(controlcard,"RATTLE").gt.0
385 c if performing umbrella sampling, fragments constrained are read from the fragment file
391 if(me.eq.king.or..not.out1file) then
393 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
395 write (iout,'(a)') "The units are:"
396 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
397 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
398 & " acceleration: angstrom/(48.9 fs)**2"
399 write (iout,'(a)') "energy: kcal/mol, temperature: K"
401 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
402 write (iout,'(a60,f10.5,a)')
403 & "Initial time step of numerical integration:",d_time,
405 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
407 write (iout,'(2a,i4,a)')
408 & "A-MTS algorithm used; initial time step for fast-varying",
409 & " short-range forces split into",ntime_split," steps."
410 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
411 & r_cut," lambda",rlamb
413 write (iout,'(2a,f10.5)')
414 & "Maximum acceleration threshold to reduce the time step",
415 & "/increase split number:",damax
416 write (iout,'(2a,f10.5)')
417 & "Maximum predicted energy drift to reduce the timestep",
418 & "/increase split number:",edriftmax
419 write (iout,'(a60,f10.5)')
420 & "Maximum velocity threshold to reduce velocities:",dvmax
421 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
422 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
423 if (rattle) write (iout,'(a60)')
424 & "Rattle algorithm used to constrain the virtual bonds"
428 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
429 call reada(controlcard,"RWAT",rwat,1.4d0)
430 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
431 surfarea=index(controlcard,"SURFAREA").gt.0
432 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
433 if(me.eq.king.or..not.out1file)then
434 write (iout,'(/a,$)') "Langevin dynamics calculation"
437 & " with direct integration of Langevin equations"
438 else if (lang.eq.2) then
439 write (iout,'(a/)') " with TINKER stochasic MD integrator"
440 else if (lang.eq.3) then
441 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
442 else if (lang.eq.4) then
443 write (iout,'(a/)') " in overdamped mode"
445 write (iout,'(//a,i5)')
446 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
449 write (iout,'(a60,f10.5)') "Temperature:",t_bath
450 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
451 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
452 write (iout,'(a60,f10.5)')
453 & "Scaling factor of the friction forces:",scal_fric
454 if (surfarea) write (iout,'(2a,i10,a)')
455 & "Friction coefficients will be scaled by solvent-accessible",
456 & " surface area every",reset_fricmat," steps."
458 c Calculate friction coefficients and bounds of stochastic forces
459 eta=6*pi*cPoise*etawat
460 if(me.eq.king.or..not.out1file)
461 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
463 gamp=scal_fric*(pstok+rwat)*eta
464 stdfp=dsqrt(2*Rb*t_bath/d_time)
466 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
467 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
469 if(me.eq.king.or..not.out1file)then
470 write (iout,'(/2a/)')
471 & "Radii of site types and friction coefficients and std's of",
472 & " stochastic forces of fully exposed sites"
473 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
475 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
476 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
480 if(me.eq.king.or..not.out1file)then
481 write (iout,'(a)') "Berendsen bath calculation"
482 write (iout,'(a60,f10.5)') "Temperature:",t_bath
483 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
485 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
486 & count_reset_moment," steps"
488 & write (iout,'(a,i10,a)')
489 & "Velocities will be reset at random every",count_reset_vel,
493 if(me.eq.king.or..not.out1file)
494 & write (iout,'(a31)') "Microcanonical mode calculation"
496 if(me.eq.king.or..not.out1file)then
497 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
499 write(iout,*) "MD running with constraints."
500 write(iout,*) "Equilibration time ", eq_time, " mtus."
501 write(iout,*) "Constraining ", nfrag," fragments."
502 write(iout,*) "Length of each fragment, weight and q0:"
504 write (iout,*) "Set of restraints #",iset
506 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
507 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
509 write(iout,*) "constraints between ", npair, "fragments."
510 write(iout,*) "constraint pairs, weights and q0:"
512 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
513 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
515 write(iout,*) "angle constraints within ", nfrag_back,
516 & "backbone fragments."
517 write(iout,*) "fragment, weights:"
519 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
520 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
521 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
524 iset=mod(kolor,nset)+1
527 if(me.eq.king.or..not.out1file)
528 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
531 c------------------------------------------------------------------------------
534 C Read molecular data.
536 implicit real*8 (a-h,o-z)
542 include 'COMMON.IOUNITS'
545 include 'COMMON.INTERACT'
546 include 'COMMON.LOCAL'
547 include 'COMMON.NAMES'
548 include 'COMMON.CHAIN'
549 include 'COMMON.FFIELD'
550 include 'COMMON.SBRIDGE'
551 include 'COMMON.HEADER'
552 include 'COMMON.CONTROL'
553 include 'COMMON.DBASE'
554 include 'COMMON.THREAD'
555 include 'COMMON.CONTACTS'
556 include 'COMMON.TORCNSTR'
557 include 'COMMON.TIME1'
558 include 'COMMON.BOUNDS'
560 include 'COMMON.SETUP'
561 character*4 sequence(maxres)
563 double precision x(maxvar)
564 character*256 pdbfile
565 character*320 weightcard
566 character*80 weightcard_t,ucase
567 dimension itype_pdb(maxres)
568 common /pizda/ itype_pdb
569 logical seq_comp,fail
570 double precision energia(0:n_ene)
576 C Read weights of the subsequent energy terms.
577 call card_concat(weightcard)
578 call reada(weightcard,'WLONG',wlong,1.0D0)
579 call reada(weightcard,'WSC',wsc,wlong)
580 call reada(weightcard,'WSCP',wscp,wlong)
581 call reada(weightcard,'WELEC',welec,1.0D0)
582 call reada(weightcard,'WVDWPP',wvdwpp,welec)
583 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
584 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
585 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
586 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
587 call reada(weightcard,'WTURN3',wturn3,1.0D0)
588 call reada(weightcard,'WTURN4',wturn4,1.0D0)
589 call reada(weightcard,'WTURN6',wturn6,1.0D0)
590 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
591 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
592 call reada(weightcard,'WBOND',wbond,1.0D0)
593 call reada(weightcard,'WTOR',wtor,1.0D0)
594 call reada(weightcard,'WTORD',wtor_d,1.0D0)
595 call reada(weightcard,'WANG',wang,1.0D0)
596 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
597 call reada(weightcard,'SCAL14',scal14,0.4D0)
598 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
599 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
600 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
601 call reada(weightcard,'TEMP0',temp0,300.0d0)
602 call reada(weightcard,'WLT',wliptran,0.0D0)
603 if (index(weightcard,'SOFT').gt.0) ipot=6
604 C 12/1/95 Added weight for the multi-body term WCORR
605 call reada(weightcard,'WCORRH',wcorr,1.0D0)
606 if (wcorr4.gt.0.0d0) wcorr=wcorr4
626 if(me.eq.king.or..not.out1file)
627 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
628 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
630 10 format (/'Energy-term weights (unscaled):'//
631 & 'WSCC= ',f10.6,' (SC-SC)'/
632 & 'WSCP= ',f10.6,' (SC-p)'/
633 & 'WELEC= ',f10.6,' (p-p electr)'/
634 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
635 & 'WBOND= ',f10.6,' (stretching)'/
636 & 'WANG= ',f10.6,' (bending)'/
637 & 'WSCLOC= ',f10.6,' (SC local)'/
638 & 'WTOR= ',f10.6,' (torsional)'/
639 & 'WTORD= ',f10.6,' (double torsional)'/
640 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
641 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
642 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
643 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
644 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
645 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
646 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
647 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
648 & 'WTURN6= ',f10.6,' (turns, 6th order)')
649 if(me.eq.king.or..not.out1file)then
650 if (wcorr4.gt.0.0d0) then
651 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
652 & 'between contact pairs of peptide groups'
653 write (iout,'(2(a,f5.3/))')
654 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
655 & 'Range of quenching the correlation terms:',2*delt_corr
656 else if (wcorr.gt.0.0d0) then
657 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
658 & 'between contact pairs of peptide groups'
660 write (iout,'(a,f8.3)')
661 & 'Scaling factor of 1,4 SC-p interactions:',scal14
662 write (iout,'(a,f8.3)')
663 & 'General scaling factor of SC-p interactions:',scalscp
665 r0_corr=cutoff_corr-delt_corr
667 aad(i,1)=scalscp*aad(i,1)
668 aad(i,2)=scalscp*aad(i,2)
669 bad(i,1)=scalscp*bad(i,1)
670 bad(i,2)=scalscp*bad(i,2)
672 call rescale_weights(t_bath)
673 if(me.eq.king.or..not.out1file)
674 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
675 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
677 22 format (/'Energy-term weights (scaled):'//
678 & 'WSCC= ',f10.6,' (SC-SC)'/
679 & 'WSCP= ',f10.6,' (SC-p)'/
680 & 'WELEC= ',f10.6,' (p-p electr)'/
681 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
682 & 'WBOND= ',f10.6,' (stretching)'/
683 & 'WANG= ',f10.6,' (bending)'/
684 & 'WSCLOC= ',f10.6,' (SC local)'/
685 & 'WTOR= ',f10.6,' (torsional)'/
686 & 'WTORD= ',f10.6,' (double torsional)'/
687 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
688 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
689 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
690 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
691 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
692 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
693 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
694 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
695 & 'WTURN6= ',f10.6,' (turns, 6th order)')
696 if(me.eq.king.or..not.out1file)
697 & write (iout,*) "Reference temperature for weights calculation:",
699 call reada(weightcard,"D0CM",d0cm,3.78d0)
700 call reada(weightcard,"AKCM",akcm,15.1d0)
701 call reada(weightcard,"AKTH",akth,11.0d0)
702 call reada(weightcard,"AKCT",akct,12.0d0)
703 call reada(weightcard,"V1SS",v1ss,-1.08d0)
704 call reada(weightcard,"V2SS",v2ss,7.61d0)
705 call reada(weightcard,"V3SS",v3ss,13.7d0)
706 call reada(weightcard,"EBR",ebr,-5.50D0)
707 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
709 dyn_ss_mask(i)=.false.
713 dyn_ssbond_ij(i,j)=1.0d300
716 call reada(weightcard,"HT",Ht,0.0D0)
718 ss_depth=ebr/wsc-0.25*eps(1,1)
719 Ht=Ht/wsc-0.25*eps(1,1)
720 akcm=akcm*wstrain/wsc
721 akth=akth*wstrain/wsc
722 akct=akct*wstrain/wsc
723 v1ss=v1ss*wstrain/wsc
724 v2ss=v2ss*wstrain/wsc
725 v3ss=v3ss*wstrain/wsc
727 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
730 if(me.eq.king.or..not.out1file) then
731 write (iout,*) "Parameters of the SS-bond potential:"
732 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
734 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
735 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
736 write (iout,*)" HT",Ht
737 print *,'indpdb=',indpdb,' pdbref=',pdbref
739 if (indpdb.gt.0 .or. pdbref) then
740 read(inp,'(a)') pdbfile
741 if(me.eq.king.or..not.out1file)
742 & write (iout,'(2a)') 'PDB data will be read from file ',
743 & pdbfile(:ilen(pdbfile))
744 open(ipdbin,file=pdbfile,status='old',err=33)
746 33 write (iout,'(a)') 'Error opening PDB file.'
749 c print *,'Begin reading pdb data'
751 c print *,'Finished reading pdb data'
752 if(me.eq.king.or..not.out1file)
753 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
754 & ' nstart_sup=',nstart_sup
756 itype_pdb(i)=itype(i)
760 nct=nstart_sup+nsup-1
761 call contact(.false.,ncont_ref,icont_ref,co)
764 if(me.eq.king.or..not.out1file)
765 & write(iout,*)'Adding sidechains'
769 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
772 do while (fail.and.nsi.le.maxsi)
773 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
776 if(fail) write(iout,*)'Adding sidechain failed for res ',
777 & i,' after ',nsi,' trials'
782 if (indpdb.eq.0) then
783 C Read sequence if not taken from the pdb file.
785 c print *,'nres=',nres
786 if (iscode.gt.0) then
787 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
789 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
791 C Convert sequence to numeric code
793 itype(i)=rescode(i,sequence(i),iscode)
795 C Assign initial virtual bond lengths
801 vbld(i+nres)=dsc(iabs(itype(i)))
802 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
803 c write (iout,*) "i",i," itype",itype(i),
804 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
808 c print '(20i4)',(itype(i),i=1,nres)
811 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
813 if (itype(i).eq.ntyp1) then
817 else if (iabs(itype(i+1)).ne.20) then
819 else if (iabs(itype(i)).ne.20) then
826 if(me.eq.king.or..not.out1file)then
827 write (iout,*) "ITEL"
829 write (iout,*) i,itype(i),itel(i)
831 print *,'Call Read_Bridge.'
834 C 8/13/98 Set limits to generating the dihedral angles
839 read (inp,*) ndih_constr
840 if (ndih_constr.gt.0) then
842 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
843 if(me.eq.king.or..not.out1file)then
845 & 'There are',ndih_constr,' constraints on phi angles.'
847 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
851 phi0(i)=deg2rad*phi0(i)
852 drange(i)=deg2rad*drange(i)
854 if(me.eq.king.or..not.out1file)
855 & write (iout,*) 'FTORS',ftors
858 phibound(1,ii) = phi0(i)-drange(i)
859 phibound(2,ii) = phi0(i)+drange(i)
866 write (iout,'(a)') 'Boundaries in phi angle sampling:'
868 write (iout,'(a3,i5,2f10.1)')
869 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
875 cd print *,'NNT=',NNT,' NCT=',NCT
876 if (itype(1).eq.ntyp1) nnt=2
877 if (itype(nres).eq.ntyp1) nct=nct-1
879 if(me.eq.king.or..not.out1file)
880 & write (iout,'(a,i3)') 'nsup=',nsup
882 if (nsup.le.(nct-nnt+1)) then
883 do i=0,nct-nnt+1-nsup
884 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
890 & 'Error - sequences to be superposed do not match.'
893 do i=0,nsup-(nct-nnt+1)
894 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
896 nstart_sup=nstart_sup+i
902 & 'Error - sequences to be superposed do not match.'
905 if (nsup.eq.0) nsup=nct-nnt
906 if (nstart_sup.eq.0) nstart_sup=nnt
907 if (nstart_seq.eq.0) nstart_seq=nnt
908 if(me.eq.king.or..not.out1file)
909 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
910 & ' nstart_seq=',nstart_seq
912 c--- Zscore rms -------
913 if (nz_start.eq.0) nz_start=nnt
914 if (nz_end.eq.0 .and. nsup.gt.0) then
916 else if (nz_end.eq.0) then
919 if(me.eq.king.or..not.out1file)then
920 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
921 write (iout,*) 'IZ_SC=',iz_sc
923 c----------------------
926 if (.not.pdbref) then
927 call read_angles(inp,*38)
929 38 write (iout,'(a)') 'Error reading reference structure.'
931 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
932 stop 'Error reading reference structure'
936 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
946 call contact(.true.,ncont_ref,icont_ref,co)
948 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
950 if (constr_dist.gt.0) call read_dist_constr
951 write (iout,*) "After read_dist_constr nhpb",nhpb
953 if(me.eq.king.or..not.out1file)
954 & write (iout,*) 'Contact order:',co
956 if(me.eq.king.or..not.out1file)
957 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
960 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
962 if(me.eq.king.or..not.out1file)
963 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
964 & icont_ref(1,i),' ',
965 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
969 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
970 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
971 & modecalc.ne.10) then
972 C If input structure hasn't been supplied from the PDB file read or generate
974 if (iranconf.eq.0 .and. .not. extconf) then
975 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
976 & write (iout,'(a)') 'Initial geometry will be read in.'
978 read(inp,'(8f10.5)',end=36,err=36)
979 & ((c(l,k),l=1,3),k=1,nres),
980 & ((c(l,k+nres),l=1,3),k=nnt,nct)
981 write (iout,*) "Exit READ_CART"
982 write (iout,'(8f10.5)')
983 & ((c(l,k),l=1,3),k=1,nres),
984 & ((c(l,k+nres),l=1,3),k=nnt,nct)
985 call int_from_cart1(.true.)
986 write (iout,*) "Finish INT_TO_CART"
989 dc(j,i)=c(j,i+1)-c(j,i)
990 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
994 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
996 dc(j,i+nres)=c(j,i+nres)-c(j,i)
997 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1003 call read_angles(inp,*36)
1006 36 write (iout,'(a)') 'Error reading angle file.'
1008 call mpi_finalize( MPI_COMM_WORLD,IERR )
1010 stop 'Error reading angle file.'
1012 else if (extconf) then
1013 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1014 & write (iout,'(a)') 'Extended chain initial geometry.'
1016 theta(i)=90d0*deg2rad
1019 phi(i)=180d0*deg2rad
1022 alph(i)=110d0*deg2rad
1025 omeg(i)=-120d0*deg2rad
1026 if (itype(i).le.0) omeg(i)=-omeg(i)
1029 if(me.eq.king.or..not.out1file)
1030 & write (iout,'(a)') 'Random-generated initial geometry.'
1034 if (me.eq.king .or. fg_rank.eq.0 .and. (
1035 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1039 call gen_rand_conf(itmp,*30)
1041 30 write (iout,*) 'Failed to generate random conformation',
1042 & ', itrial=',itrial
1043 write (*,*) 'Processor:',me,
1044 & ' Failed to generate random conformation',
1054 write (iout,'(a,i3,a)') 'Processor:',me,
1055 & ' error in generating random conformation.'
1056 write (*,'(a,i3,a)') 'Processor:',me,
1057 & ' error in generating random conformation.'
1060 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1065 & ' error in generating random conformation.'
1070 elseif (modecalc.eq.4) then
1071 read (inp,'(a)') intinname
1072 open (intin,file=intinname,status='old',err=333)
1073 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1074 & write (iout,'(a)') 'intinname',intinname
1075 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1077 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1079 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1081 stop 'Error opening angle file.'
1085 C Generate distance constraints, if the PDB structure is to be regularized.
1086 if (nthread.gt.0) then
1087 call read_threadbase
1090 if (me.eq.king .or. .not. out1file)
1092 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1093 write (iout,'(/a,i3,a)')
1094 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1095 write (iout,'(20i4)') (iss(i),i=1,ns)
1097 write(iout,*)"Running with dynamic disulfide-bond formation"
1099 write (iout,'(/a/)') 'Pre-formed links are:'
1105 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1106 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1112 if (ns.gt.0.and.dyn_ss) then
1116 forcon(i-nss)=forcon(i)
1123 dyn_ss_mask(iss(i))=.true.
1126 if (i2ndstr.gt.0) call secstrp2dihc
1127 c call geom_to_var(nvar,x)
1128 c call etotal(energia(0))
1129 c call enerprint(energia(0))
1130 c call briefout(0,etot)
1132 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1133 cd write (iout,'(a)') 'Variable list:'
1134 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1136 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1137 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1138 & 'Processor',myrank,': end reading molecular data.'
1142 c--------------------------------------------------------------------------
1143 logical function seq_comp(itypea,itypeb,length)
1145 integer length,itypea(length),itypeb(length)
1148 if (itypea(i).ne.itypeb(i)) then
1156 c-----------------------------------------------------------------------------
1157 subroutine read_bridge
1158 C Read information about disulfide bridges.
1159 implicit real*8 (a-h,o-z)
1160 include 'DIMENSIONS'
1164 include 'COMMON.IOUNITS'
1165 include 'COMMON.GEO'
1166 include 'COMMON.VAR'
1167 include 'COMMON.INTERACT'
1168 include 'COMMON.LOCAL'
1169 include 'COMMON.NAMES'
1170 include 'COMMON.CHAIN'
1171 include 'COMMON.FFIELD'
1172 include 'COMMON.SBRIDGE'
1173 include 'COMMON.HEADER'
1174 include 'COMMON.CONTROL'
1175 include 'COMMON.DBASE'
1176 include 'COMMON.THREAD'
1177 include 'COMMON.TIME1'
1178 include 'COMMON.SETUP'
1179 C Read bridging residues.
1180 read (inp,*) ns,(iss(i),i=1,ns)
1182 if(me.eq.king.or..not.out1file)
1183 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1184 C Check whether the specified bridging residues are cystines.
1186 if (itype(iss(i)).ne.1) then
1187 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1188 & 'Do you REALLY think that the residue ',
1189 & restyp(itype(iss(i))),i,
1190 & ' can form a disulfide bridge?!!!'
1191 write (*,'(2a,i3,a)')
1192 & 'Do you REALLY think that the residue ',
1193 & restyp(itype(iss(i))),i,
1194 & ' can form a disulfide bridge?!!!'
1196 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1201 C Read preformed bridges.
1203 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1205 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1208 C Check if the residues involved in bridges are in the specified list of
1209 C bridging residues.
1212 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1213 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1214 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1215 & ' contains residues present in other pairs.'
1216 write (*,'(a,i3,a)') 'Disulfide pair',i,
1217 & ' contains residues present in other pairs.'
1219 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1225 if (ihpb(i).eq.iss(j)) goto 10
1227 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1230 if (jhpb(i).eq.iss(j)) goto 20
1232 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1238 ihpb(i)=ihpb(i)+nres
1239 jhpb(i)=jhpb(i)+nres
1245 c----------------------------------------------------------------------------
1246 subroutine read_x(kanal,*)
1247 implicit real*8 (a-h,o-z)
1248 include 'DIMENSIONS'
1249 include 'COMMON.GEO'
1250 include 'COMMON.VAR'
1251 include 'COMMON.CHAIN'
1252 include 'COMMON.IOUNITS'
1253 include 'COMMON.CONTROL'
1254 include 'COMMON.LOCAL'
1255 include 'COMMON.INTERACT'
1256 c Read coordinates from input
1258 read(kanal,'(8f10.5)',end=10,err=10)
1259 & ((c(l,k),l=1,3),k=1,nres),
1260 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1263 c(j,2*nres)=c(j,nres)
1265 call int_from_cart1(.false.)
1268 dc(j,i)=c(j,i+1)-c(j,i)
1269 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1273 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1275 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1276 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1284 c----------------------------------------------------------------------------
1285 subroutine read_threadbase
1286 implicit real*8 (a-h,o-z)
1287 include 'DIMENSIONS'
1288 include 'COMMON.IOUNITS'
1289 include 'COMMON.GEO'
1290 include 'COMMON.VAR'
1291 include 'COMMON.INTERACT'
1292 include 'COMMON.LOCAL'
1293 include 'COMMON.NAMES'
1294 include 'COMMON.CHAIN'
1295 include 'COMMON.FFIELD'
1296 include 'COMMON.SBRIDGE'
1297 include 'COMMON.HEADER'
1298 include 'COMMON.CONTROL'
1299 include 'COMMON.DBASE'
1300 include 'COMMON.THREAD'
1301 include 'COMMON.TIME1'
1302 C Read pattern database for threading.
1303 read (icbase,*) nseq
1305 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1306 & nres_base(2,i),nres_base(3,i)
1307 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1309 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1310 c & nres_base(2,i),nres_base(3,i)
1311 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1315 if (weidis.eq.0.0D0) weidis=0.1D0
1324 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1325 write (iout,'(a,i5)') 'nexcl: ',nexcl
1326 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1329 c------------------------------------------------------------------------------
1330 subroutine setup_var
1331 implicit real*8 (a-h,o-z)
1332 include 'DIMENSIONS'
1333 include 'COMMON.IOUNITS'
1334 include 'COMMON.GEO'
1335 include 'COMMON.VAR'
1336 include 'COMMON.INTERACT'
1337 include 'COMMON.LOCAL'
1338 include 'COMMON.NAMES'
1339 include 'COMMON.CHAIN'
1340 include 'COMMON.FFIELD'
1341 include 'COMMON.SBRIDGE'
1342 include 'COMMON.HEADER'
1343 include 'COMMON.CONTROL'
1344 include 'COMMON.DBASE'
1345 include 'COMMON.THREAD'
1346 include 'COMMON.TIME1'
1347 C Set up variable list.
1353 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1355 ialph(i,1)=nvar+nside
1359 if (indphi.gt.0) then
1361 else if (indback.gt.0) then
1366 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1369 c----------------------------------------------------------------------------
1370 subroutine gen_dist_constr
1371 C Generate CA distance constraints.
1372 implicit real*8 (a-h,o-z)
1373 include 'DIMENSIONS'
1374 include 'COMMON.IOUNITS'
1375 include 'COMMON.GEO'
1376 include 'COMMON.VAR'
1377 include 'COMMON.INTERACT'
1378 include 'COMMON.LOCAL'
1379 include 'COMMON.NAMES'
1380 include 'COMMON.CHAIN'
1381 include 'COMMON.FFIELD'
1382 include 'COMMON.SBRIDGE'
1383 include 'COMMON.HEADER'
1384 include 'COMMON.CONTROL'
1385 include 'COMMON.DBASE'
1386 include 'COMMON.THREAD'
1387 include 'COMMON.TIME1'
1388 dimension itype_pdb(maxres)
1389 common /pizda/ itype_pdb
1391 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1392 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1393 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1395 do i=nstart_sup,nstart_sup+nsup-1
1396 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1397 cd & ' seq_pdb', restyp(itype_pdb(i))
1398 do j=i+2,nstart_sup+nsup-1
1400 ihpb(nhpb)=i+nstart_seq-nstart_sup
1401 jhpb(nhpb)=j+nstart_seq-nstart_sup
1403 dhpb(nhpb)=dist(i,j)
1406 cd write (iout,'(a)') 'Distance constraints:'
1411 cd if (ii.gt.nres) then
1416 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1417 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1418 cd & dhpb(i),forcon(i)
1422 c----------------------------------------------------------------------------
1424 implicit real*8 (a-h,o-z)
1425 include 'DIMENSIONS'
1426 include 'COMMON.MAP'
1427 include 'COMMON.IOUNITS'
1428 character*3 angid(4) /'THE','PHI','ALP','OME'/
1429 character*80 mapcard,ucase
1431 read (inp,'(a)') mapcard
1432 mapcard=ucase(mapcard)
1433 if (index(mapcard,'PHI').gt.0) then
1435 else if (index(mapcard,'THE').gt.0) then
1437 else if (index(mapcard,'ALP').gt.0) then
1439 else if (index(mapcard,'OME').gt.0) then
1442 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1443 stop 'Error - illegal variable spec in MAP card.'
1445 call readi (mapcard,'RES1',res1(imap),0)
1446 call readi (mapcard,'RES2',res2(imap),0)
1447 if (res1(imap).eq.0) then
1448 res1(imap)=res2(imap)
1449 else if (res2(imap).eq.0) then
1450 res2(imap)=res1(imap)
1452 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1454 & 'Error - illegal definition of variable group in MAP.'
1455 stop 'Error - illegal definition of variable group in MAP.'
1457 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1458 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1459 call readi(mapcard,'NSTEP',nstep(imap),0)
1460 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1462 & 'Illegal boundary and/or step size specification in MAP.'
1463 stop 'Illegal boundary and/or step size specification in MAP.'
1468 c----------------------------------------------------------------------------
1470 implicit real*8 (a-h,o-z)
1471 include 'DIMENSIONS'
1472 include 'COMMON.IOUNITS'
1473 include 'COMMON.GEO'
1474 include 'COMMON.CSA'
1475 include 'COMMON.BANK'
1476 include 'COMMON.CONTROL'
1478 character*620 mcmcard
1479 call card_concat(mcmcard)
1481 call readi(mcmcard,'NCONF',nconf,50)
1482 call readi(mcmcard,'NADD',nadd,0)
1483 call readi(mcmcard,'JSTART',jstart,1)
1484 call readi(mcmcard,'JEND',jend,1)
1485 call readi(mcmcard,'NSTMAX',nstmax,500000)
1486 call readi(mcmcard,'N0',n0,1)
1487 call readi(mcmcard,'N1',n1,6)
1488 call readi(mcmcard,'N2',n2,4)
1489 call readi(mcmcard,'N3',n3,0)
1490 call readi(mcmcard,'N4',n4,0)
1491 call readi(mcmcard,'N5',n5,0)
1492 call readi(mcmcard,'N6',n6,10)
1493 call readi(mcmcard,'N7',n7,0)
1494 call readi(mcmcard,'N8',n8,0)
1495 call readi(mcmcard,'N9',n9,0)
1496 call readi(mcmcard,'N14',n14,0)
1497 call readi(mcmcard,'N15',n15,0)
1498 call readi(mcmcard,'N16',n16,0)
1499 call readi(mcmcard,'N17',n17,0)
1500 call readi(mcmcard,'N18',n18,0)
1502 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1504 call readi(mcmcard,'NDIFF',ndiff,2)
1505 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1506 call readi(mcmcard,'IS1',is1,1)
1507 call readi(mcmcard,'IS2',is2,8)
1508 call readi(mcmcard,'NRAN0',nran0,4)
1509 call readi(mcmcard,'NRAN1',nran1,2)
1510 call readi(mcmcard,'IRR',irr,1)
1511 call readi(mcmcard,'NSEED',nseed,20)
1512 call readi(mcmcard,'NTOTAL',ntotal,10000)
1513 call reada(mcmcard,'CUT1',cut1,2.0d0)
1514 call reada(mcmcard,'CUT2',cut2,5.0d0)
1515 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1516 call readi(mcmcard,'ICMAX',icmax,3)
1517 call readi(mcmcard,'IRESTART',irestart,0)
1518 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1521 call reada(mcmcard,'DELE',dele,20.0d0)
1522 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1523 call readi(mcmcard,'IREF',iref,0)
1524 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1525 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1526 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1527 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1528 write (iout,*) "NCONF_IN",nconf_in
1531 c----------------------------------------------------------------------------
1532 cfmc subroutine mcmfread
1533 cfmc implicit real*8 (a-h,o-z)
1534 cfmc include 'DIMENSIONS'
1535 cfmc include 'COMMON.MCMF'
1536 cfmc include 'COMMON.IOUNITS'
1537 cfmc include 'COMMON.GEO'
1538 cfmc character*80 ucase
1539 cfmc character*620 mcmcard
1540 cfmc call card_concat(mcmcard)
1542 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1543 cfmc write(iout,*)'MAXRANT=',maxrant
1544 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1545 cfmc write(iout,*)'MAXFAM=',maxfam
1546 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1547 cfmc write(iout,*)'NNET1=',nnet1
1548 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1549 cfmc write(iout,*)'NNET2=',nnet2
1550 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1551 cfmc write(iout,*)'NNET3=',nnet3
1552 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1553 cfmc write(iout,*)'ILASTT=',ilastt
1554 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1555 cfmc write(iout,*)'MAXSTR=',maxstr
1556 cfmc maxstr_f=maxstr/maxfam
1557 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1558 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1559 cfmc write(iout,*)'NMCMF=',nmcmf
1560 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1561 cfmc write(iout,*)'IFOCUS=',ifocus
1562 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1563 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1564 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1565 cfmc write(iout,*)'INTPRT=',intprt
1566 cfmc call readi(mcmcard,'IPRT',iprt,100)
1567 cfmc write(iout,*)'IPRT=',iprt
1568 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1569 cfmc write(iout,*)'IMAXTR=',imaxtr
1570 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1571 cfmc write(iout,*)'MAXEVEN=',maxeven
1572 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1573 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1574 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1575 cfmc write(iout,*)'INIMIN=',inimin
1576 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1577 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1578 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1579 cfmc write(iout,*)'NTHREAD=',nthread
1580 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1581 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1582 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1583 cfmc write(iout,*)'MAXPERT=',maxpert
1584 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1585 cfmc write(iout,*)'IRMSD=',irmsd
1586 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1587 cfmc write(iout,*)'DENEMIN=',denemin
1588 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1589 cfmc write(iout,*)'RCUT1S=',rcut1s
1590 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1591 cfmc write(iout,*)'RCUT1E=',rcut1e
1592 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1593 cfmc write(iout,*)'RCUT2S=',rcut2s
1594 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1595 cfmc write(iout,*)'RCUT2E=',rcut2e
1596 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1597 cfmc write(iout,*)'DPERT1=',d_pert1
1598 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1599 cfmc write(iout,*)'DPERT1A=',d_pert1a
1600 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1601 cfmc write(iout,*)'DPERT2=',d_pert2
1602 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1603 cfmc write(iout,*)'DPERT2A=',d_pert2a
1604 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1605 cfmc write(iout,*)'DPERT2B=',d_pert2b
1606 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1607 cfmc write(iout,*)'DPERT2C=',d_pert2c
1608 cfmc d_pert1=deg2rad*d_pert1
1609 cfmc d_pert1a=deg2rad*d_pert1a
1610 cfmc d_pert2=deg2rad*d_pert2
1611 cfmc d_pert2a=deg2rad*d_pert2a
1612 cfmc d_pert2b=deg2rad*d_pert2b
1613 cfmc d_pert2c=deg2rad*d_pert2c
1614 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1615 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1616 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1617 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1618 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1619 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1620 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1621 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1622 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1623 cfmc write(iout,*)'RCUTINI=',rcutini
1624 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1625 cfmc write(iout,*)'GRAT=',grat
1626 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1627 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1631 c----------------------------------------------------------------------------
1633 implicit real*8 (a-h,o-z)
1634 include 'DIMENSIONS'
1635 include 'COMMON.MCM'
1636 include 'COMMON.MCE'
1637 include 'COMMON.IOUNITS'
1639 character*320 mcmcard
1640 call card_concat(mcmcard)
1641 call readi(mcmcard,'MAXACC',maxacc,100)
1642 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1643 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1644 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1645 call readi(mcmcard,'MAXREPM',maxrepm,200)
1646 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1647 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1648 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1649 call reada(mcmcard,'E_UP',e_up,5.0D0)
1650 call reada(mcmcard,'DELTE',delte,0.1D0)
1651 call readi(mcmcard,'NSWEEP',nsweep,5)
1652 call readi(mcmcard,'NSTEPH',nsteph,0)
1653 call readi(mcmcard,'NSTEPC',nstepc,0)
1654 call reada(mcmcard,'TMIN',tmin,298.0D0)
1655 call reada(mcmcard,'TMAX',tmax,298.0D0)
1656 call readi(mcmcard,'NWINDOW',nwindow,0)
1657 call readi(mcmcard,'PRINT_MC',print_mc,0)
1658 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1659 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1660 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1661 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1662 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1663 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1664 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1665 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1666 if (nwindow.gt.0) then
1667 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1669 winlen(i)=winend(i)-winstart(i)+1
1672 if (tmax.lt.tmin) tmax=tmin
1673 if (tmax.eq.tmin) then
1677 if (nstepc.gt.0 .and. nsteph.gt.0) then
1678 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1679 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1681 C Probabilities of different move types
1682 sumpro_type(0)=0.0D0
1683 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1684 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1685 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1686 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1687 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1688 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1689 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1691 print *,'i',i,' sumprotype',sumpro_type(i)
1692 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1693 print *,'i',i,' sumprotype',sumpro_type(i)
1697 c----------------------------------------------------------------------------
1698 subroutine read_minim
1699 implicit real*8 (a-h,o-z)
1700 include 'DIMENSIONS'
1701 include 'COMMON.MINIM'
1702 include 'COMMON.IOUNITS'
1704 character*320 minimcard
1705 call card_concat(minimcard)
1706 call readi(minimcard,'MAXMIN',maxmin,2000)
1707 call readi(minimcard,'MAXFUN',maxfun,5000)
1708 call readi(minimcard,'MINMIN',minmin,maxmin)
1709 call readi(minimcard,'MINFUN',minfun,maxmin)
1710 call reada(minimcard,'TOLF',tolf,1.0D-2)
1711 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1712 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1713 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1714 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1715 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1716 & 'Options in energy minimization:'
1717 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1718 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1719 & 'MinMin:',MinMin,' MinFun:',MinFun,
1720 & ' TolF:',TolF,' RTolF:',RTolF
1723 c----------------------------------------------------------------------------
1724 subroutine read_angles(kanal,*)
1725 implicit real*8 (a-h,o-z)
1726 include 'DIMENSIONS'
1727 include 'COMMON.GEO'
1728 include 'COMMON.VAR'
1729 include 'COMMON.CHAIN'
1730 include 'COMMON.IOUNITS'
1731 include 'COMMON.CONTROL'
1732 c Read angles from input
1734 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1735 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1736 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1737 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1740 c 9/7/01 avoid 180 deg valence angle
1741 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1743 theta(i)=deg2rad*theta(i)
1744 phi(i)=deg2rad*phi(i)
1745 alph(i)=deg2rad*alph(i)
1746 omeg(i)=deg2rad*omeg(i)
1751 c----------------------------------------------------------------------------
1752 subroutine reada(rekord,lancuch,wartosc,default)
1754 character*(*) rekord,lancuch
1755 double precision wartosc,default
1758 iread=index(rekord,lancuch)
1759 if (iread.eq.0) then
1763 iread=iread+ilen(lancuch)+1
1764 read (rekord(iread:),*,err=10,end=10) wartosc
1769 c----------------------------------------------------------------------------
1770 subroutine readi(rekord,lancuch,wartosc,default)
1772 character*(*) rekord,lancuch
1773 integer wartosc,default
1776 iread=index(rekord,lancuch)
1777 if (iread.eq.0) then
1781 iread=iread+ilen(lancuch)+1
1782 read (rekord(iread:),*,err=10,end=10) wartosc
1787 c----------------------------------------------------------------------------
1788 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1791 integer tablica(dim),default
1792 character*(*) rekord,lancuch
1799 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1800 if (iread.eq.0) return
1801 iread=iread+ilen(lancuch)+1
1802 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1805 c----------------------------------------------------------------------------
1806 subroutine multreada(rekord,lancuch,tablica,dim,default)
1809 double precision tablica(dim),default
1810 character*(*) rekord,lancuch
1817 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1818 if (iread.eq.0) return
1819 iread=iread+ilen(lancuch)+1
1820 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1823 c----------------------------------------------------------------------------
1824 subroutine openunits
1825 implicit real*8 (a-h,o-z)
1826 include 'DIMENSIONS'
1829 character*16 form,nodename
1832 include 'COMMON.SETUP'
1833 include 'COMMON.IOUNITS'
1835 include 'COMMON.CONTROL'
1836 integer lenpre,lenpot,ilen,lentmp
1838 character*3 out1file_text,ucase
1841 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1842 call getenv_loc("PREFIX",prefix)
1844 call getenv_loc("POT",pot)
1845 call getenv_loc("DIRTMP",tmpdir)
1846 call getenv_loc("CURDIR",curdir)
1847 call getenv_loc("OUT1FILE",out1file_text)
1848 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1849 out1file_text=ucase(out1file_text)
1850 if (out1file_text(1:1).eq."Y") then
1853 out1file=fg_rank.gt.0
1858 if (lentmp.gt.0) then
1859 write (*,'(80(1h!))')
1860 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1861 write (*,'(80(1h!))')
1862 write (*,*)"All output files will be on node /tmp directory."
1864 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1865 if (me.eq.king) then
1866 write (*,*) "The master node is ",nodename
1867 else if (fg_rank.eq.0) then
1868 write (*,*) "I am the CG slave node ",nodename
1870 write (*,*) "I am the FG slave node ",nodename
1873 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1874 lenpre = lentmp+lenpre+1
1876 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1877 C Get the names and open the input files
1878 #if defined(WINIFL) || defined(WINPGI)
1879 open(1,file=pref_orig(:ilen(pref_orig))//
1880 & '.inp',status='old',readonly,shared)
1881 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1882 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1883 C Get parameter filenames and open the parameter files.
1884 call getenv_loc('BONDPAR',bondname)
1885 open (ibond,file=bondname,status='old',readonly,shared)
1886 call getenv_loc('THETPAR',thetname)
1887 open (ithep,file=thetname,status='old',readonly,shared)
1888 call getenv_loc('ROTPAR',rotname)
1889 open (irotam,file=rotname,status='old',readonly,shared)
1890 call getenv_loc('TORPAR',torname)
1891 open (itorp,file=torname,status='old',readonly,shared)
1892 call getenv_loc('TORDPAR',tordname)
1893 open (itordp,file=tordname,status='old',readonly,shared)
1894 call getenv_loc('FOURIER',fouriername)
1895 open (ifourier,file=fouriername,status='old',readonly,shared)
1896 call getenv_loc('ELEPAR',elename)
1897 open (ielep,file=elename,status='old',readonly,shared)
1898 call getenv_loc('SIDEPAR',sidename)
1899 open (isidep,file=sidename,status='old',readonly,shared)
1900 #elif (defined CRAY) || (defined AIX)
1901 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1903 c print *,"Processor",myrank," opened file 1"
1904 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1905 c print *,"Processor",myrank," opened file 9"
1906 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1907 C Get parameter filenames and open the parameter files.
1908 call getenv_loc('BONDPAR',bondname)
1909 open (ibond,file=bondname,status='old',action='read')
1910 c print *,"Processor",myrank," opened file IBOND"
1911 call getenv_loc('THETPAR',thetname)
1912 open (ithep,file=thetname,status='old',action='read')
1913 c print *,"Processor",myrank," opened file ITHEP"
1914 call getenv_loc('ROTPAR',rotname)
1915 open (irotam,file=rotname,status='old',action='read')
1916 c print *,"Processor",myrank," opened file IROTAM"
1917 call getenv_loc('TORPAR',torname)
1918 open (itorp,file=torname,status='old',action='read')
1919 c print *,"Processor",myrank," opened file ITORP"
1920 call getenv_loc('TORDPAR',tordname)
1921 open (itordp,file=tordname,status='old',action='read')
1922 c print *,"Processor",myrank," opened file ITORDP"
1923 call getenv_loc('SCCORPAR',sccorname)
1924 open (isccor,file=sccorname,status='old',action='read')
1925 c print *,"Processor",myrank," opened file ISCCOR"
1926 call getenv_loc('FOURIER',fouriername)
1927 open (ifourier,file=fouriername,status='old',action='read')
1928 c print *,"Processor",myrank," opened file IFOURIER"
1929 call getenv_loc('ELEPAR',elename)
1930 open (ielep,file=elename,status='old',action='read')
1931 c print *,"Processor",myrank," opened file IELEP"
1932 call getenv_loc('SIDEPAR',sidename)
1933 open (isidep,file=sidename,status='old',action='read')
1934 c print *,"Processor",myrank," opened file ISIDEP"
1935 c print *,"Processor",myrank," opened parameter files"
1937 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1938 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1939 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1940 C Get parameter filenames and open the parameter files.
1941 call getenv_loc('BONDPAR',bondname)
1942 open (ibond,file=bondname,status='old')
1943 call getenv_loc('THETPAR',thetname)
1944 open (ithep,file=thetname,status='old')
1945 call getenv_loc('ROTPAR',rotname)
1946 open (irotam,file=rotname,status='old')
1947 call getenv_loc('TORPAR',torname)
1948 open (itorp,file=torname,status='old')
1949 call getenv_loc('TORDPAR',tordname)
1950 open (itordp,file=tordname,status='old')
1951 call getenv_loc('SCCORPAR',sccorname)
1952 open (isccor,file=sccorname,status='old')
1953 call getenv_loc('FOURIER',fouriername)
1954 open (ifourier,file=fouriername,status='old')
1955 call getenv_loc('ELEPAR',elename)
1956 open (ielep,file=elename,status='old')
1957 call getenv_loc('SIDEPAR',sidename)
1958 open (isidep,file=sidename,status='old')
1960 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1962 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1963 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1964 C Get parameter filenames and open the parameter files.
1965 call getenv_loc('BONDPAR',bondname)
1966 open (ibond,file=bondname,status='old',readonly)
1967 call getenv_loc('THETPAR',thetname)
1968 open (ithep,file=thetname,status='old',readonly)
1969 call getenv_loc('ROTPAR',rotname)
1970 open (irotam,file=rotname,status='old',readonly)
1971 call getenv_loc('TORPAR',torname)
1972 open (itorp,file=torname,status='old',readonly)
1973 call getenv_loc('TORDPAR',tordname)
1974 open (itordp,file=tordname,status='old',readonly)
1975 call getenv_loc('SCCORPAR',sccorname)
1976 open (isccor,file=sccorname,status='old',readonly)
1978 call getenv_loc('THETPARPDB',thetname_pdb)
1979 print *,"thetname_pdb ",thetname_pdb
1980 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1981 print *,ithep_pdb," opened"
1983 call getenv_loc('FOURIER',fouriername)
1984 open (ifourier,file=fouriername,status='old',readonly)
1985 call getenv_loc('ELEPAR',elename)
1986 open (ielep,file=elename,status='old',readonly)
1987 call getenv_loc('SIDEPAR',sidename)
1988 open (isidep,file=sidename,status='old',readonly)
1989 call getenv_loc('LIPTRANPAR',liptranname)
1990 open (iliptranpar,file=liptranname,status='old',action='read')
1992 call getenv_loc('ROTPARPDB',rotname_pdb)
1993 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1998 C 8/9/01 In the newest version SCp interaction constants are read from a file
1999 C Use -DOLDSCP to use hard-coded constants instead.
2001 call getenv_loc('SCPPAR',scpname)
2002 #if defined(WINIFL) || defined(WINPGI)
2003 open (iscpp,file=scpname,status='old',readonly,shared)
2004 #elif (defined CRAY) || (defined AIX)
2005 open (iscpp,file=scpname,status='old',action='read')
2007 open (iscpp,file=scpname,status='old')
2009 open (iscpp,file=scpname,status='old',readonly)
2012 call getenv_loc('PATTERN',patname)
2013 #if defined(WINIFL) || defined(WINPGI)
2014 open (icbase,file=patname,status='old',readonly,shared)
2015 #elif (defined CRAY) || (defined AIX)
2016 open (icbase,file=patname,status='old',action='read')
2018 open (icbase,file=patname,status='old')
2020 open (icbase,file=patname,status='old',readonly)
2023 C Open output file only for CG processes
2024 c print *,"Processor",myrank," fg_rank",fg_rank
2025 if (fg_rank.eq.0) then
2027 if (nodes.eq.1) then
2030 npos = dlog10(dfloat(nodes-1))+1
2032 if (npos.lt.3) npos=3
2033 write (liczba,'(i1)') npos
2034 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2036 write (liczba,form) me
2037 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2038 & liczba(:ilen(liczba))
2039 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2041 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2043 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2044 & liczba(:ilen(liczba))//'.mol2'
2045 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2046 & liczba(:ilen(liczba))//'.stat'
2048 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2049 & //liczba(:ilen(liczba))//'.stat')
2050 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2053 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2054 & liczba(:ilen(liczba))//'.const'
2059 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2060 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2061 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2062 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2063 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2065 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2067 rest2name=prefix(:ilen(prefix))//'.rst'
2069 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2072 #if defined(AIX) || defined(PGI)
2073 if (me.eq.king .or. .not. out1file)
2074 & open(iout,file=outname,status='unknown')
2076 if (fg_rank.gt.0) then
2077 write (liczba,'(i3.3)') myrank/nfgtasks
2078 write (ll,'(bz,i3.3)') fg_rank
2079 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2084 open(igeom,file=intname,status='unknown',position='append')
2085 open(ipdb,file=pdbname,status='unknown')
2086 open(imol2,file=mol2name,status='unknown')
2087 open(istat,file=statname,status='unknown',position='append')
2089 c1out open(iout,file=outname,status='unknown')
2092 if (me.eq.king .or. .not.out1file)
2093 & open(iout,file=outname,status='unknown')
2095 if (fg_rank.gt.0) then
2096 write (liczba,'(i3.3)') myrank/nfgtasks
2097 write (ll,'(bz,i3.3)') fg_rank
2098 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2103 open(igeom,file=intname,status='unknown',access='append')
2104 open(ipdb,file=pdbname,status='unknown')
2105 open(imol2,file=mol2name,status='unknown')
2106 open(istat,file=statname,status='unknown',access='append')
2108 c1out open(iout,file=outname,status='unknown')
2111 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2112 csa_seed=prefix(:lenpre)//'.CSA.seed'
2113 csa_history=prefix(:lenpre)//'.CSA.history'
2114 csa_bank=prefix(:lenpre)//'.CSA.bank'
2115 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2116 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2117 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2118 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2119 csa_int=prefix(:lenpre)//'.int'
2120 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2121 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2122 csa_in=prefix(:lenpre)//'.CSA.in'
2123 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2126 write (iout,'(80(1h-))')
2127 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2128 write (iout,'(80(1h-))')
2129 write (iout,*) "Input file : ",
2130 & pref_orig(:ilen(pref_orig))//'.inp'
2131 write (iout,*) "Output file : ",
2132 & outname(:ilen(outname))
2134 write (iout,*) "Sidechain potential file : ",
2135 & sidename(:ilen(sidename))
2137 write (iout,*) "SCp potential file : ",
2138 & scpname(:ilen(scpname))
2140 write (iout,*) "Electrostatic potential file : ",
2141 & elename(:ilen(elename))
2142 write (iout,*) "Cumulant coefficient file : ",
2143 & fouriername(:ilen(fouriername))
2144 write (iout,*) "Torsional parameter file : ",
2145 & torname(:ilen(torname))
2146 write (iout,*) "Double torsional parameter file : ",
2147 & tordname(:ilen(tordname))
2148 write (iout,*) "SCCOR parameter file : ",
2149 & sccorname(:ilen(sccorname))
2150 write (iout,*) "Bond & inertia constant file : ",
2151 & bondname(:ilen(bondname))
2152 write (iout,*) "Bending parameter file : ",
2153 & thetname(:ilen(thetname))
2154 write (iout,*) "Rotamer parameter file : ",
2155 & rotname(:ilen(rotname))
2156 write (iout,*) "Thetpdb parameter file : ",
2157 & thetname_pdb(:ilen(thetname_pdb))
2158 write (iout,*) "Threading database : ",
2159 & patname(:ilen(patname))
2161 &write (iout,*)" DIRTMP : ",
2163 write (iout,'(80(1h-))')
2167 c----------------------------------------------------------------------------
2168 subroutine card_concat(card)
2169 implicit real*8 (a-h,o-z)
2170 include 'DIMENSIONS'
2171 include 'COMMON.IOUNITS'
2173 character*80 karta,ucase
2175 read (inp,'(a)') karta
2178 do while (karta(80:80).eq.'&')
2179 card=card(:ilen(card)+1)//karta(:79)
2180 read (inp,'(a)') karta
2183 card=card(:ilen(card)+1)//karta
2186 c----------------------------------------------------------------------------------
2188 implicit real*8 (a-h,o-z)
2189 include 'DIMENSIONS'
2190 include 'COMMON.CHAIN'
2191 include 'COMMON.IOUNITS'
2193 open(irest2,file=rest2name,status='unknown')
2194 read(irest2,*) totT,EK,potE,totE,t_bath
2196 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2199 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2202 read (irest2,*) iset
2207 c---------------------------------------------------------------------------------
2208 subroutine read_fragments
2209 implicit real*8 (a-h,o-z)
2210 include 'DIMENSIONS'
2214 include 'COMMON.SETUP'
2215 include 'COMMON.CHAIN'
2216 include 'COMMON.IOUNITS'
2218 include 'COMMON.CONTROL'
2219 read(inp,*) nset,nfrag,npair,nfrag_back
2220 if(me.eq.king.or..not.out1file)
2221 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2222 & " nfrag_back",nfrag_back
2224 read(inp,*) mset(iset)
2226 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2228 if(me.eq.king.or..not.out1file)
2229 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2230 & ifrag(2,i,iset), qinfrag(i,iset)
2233 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2235 if(me.eq.king.or..not.out1file)
2236 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2237 & ipair(2,i,iset), qinpair(i,iset)
2240 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2241 & wfrag_back(3,i,iset),
2242 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2243 if(me.eq.king.or..not.out1file)
2244 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2245 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2250 c-------------------------------------------------------------------------------
2251 subroutine read_dist_constr
2252 implicit real*8 (a-h,o-z)
2253 include 'DIMENSIONS'
2257 include 'COMMON.SETUP'
2258 include 'COMMON.CONTROL'
2259 include 'COMMON.CHAIN'
2260 include 'COMMON.IOUNITS'
2261 include 'COMMON.SBRIDGE'
2262 integer ifrag_(2,100),ipair_(2,100)
2263 double precision wfrag_(100),wpair_(100)
2264 character*500 controlcard
2265 c write (iout,*) "Calling read_dist_constr"
2266 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2268 call card_concat(controlcard)
2269 call readi(controlcard,"NFRAG",nfrag_,0)
2270 call readi(controlcard,"NPAIR",npair_,0)
2271 call readi(controlcard,"NDIST",ndist_,0)
2272 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2273 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2274 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2275 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2276 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2277 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2278 c write (iout,*) "IFRAG"
2280 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2282 c write (iout,*) "IPAIR"
2284 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2288 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2289 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2290 & ifrag_(2,i)=nstart_sup+nsup-1
2291 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2293 if (wfrag_(i).gt.0.0d0) then
2294 do j=ifrag_(1,i),ifrag_(2,i)-1
2295 do k=j+1,ifrag_(2,i)
2296 c write (iout,*) "j",j," k",k
2298 if (constr_dist.eq.1) then
2303 forcon(nhpb)=wfrag_(i)
2304 else if (constr_dist.eq.2) then
2305 if (ddjk.le.dist_cut) then
2310 forcon(nhpb)=wfrag_(i)
2317 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2320 if (.not.out1file .or. me.eq.king)
2321 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2322 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2324 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2325 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2332 if (wpair_(i).gt.0.0d0) then
2340 do j=ifrag_(1,ii),ifrag_(2,ii)
2341 do k=ifrag_(1,jj),ifrag_(2,jj)
2345 forcon(nhpb)=wpair_(i)
2346 dhpb(nhpb)=dist(j,k)
2348 if (.not.out1file .or. me.eq.king)
2349 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2350 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2352 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2353 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2360 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2361 if (forcon(nhpb+1).gt.0.0d0) then
2363 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2365 if (.not.out1file .or. me.eq.king)
2366 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2367 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2369 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2370 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2377 c-------------------------------------------------------------------------------
2379 subroutine flush(iu)
2384 subroutine flush(iu)
2389 c------------------------------------------------------------------------------
2390 subroutine copy_to_tmp(source)
2391 include "DIMENSIONS"
2392 include "COMMON.IOUNITS"
2393 character*(*) source
2394 character* 256 tmpfile
2398 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2399 inquire(file=tmpfile,exist=ex)
2401 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2402 & " to temporary directory..."
2403 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2404 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2408 c------------------------------------------------------------------------------
2409 subroutine move_from_tmp(source)
2410 include "DIMENSIONS"
2411 include "COMMON.IOUNITS"
2412 character*(*) source
2415 write (*,*) "Moving ",source(:ilen(source)),
2416 & " from temporary directory to working directory"
2417 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2418 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2421 c------------------------------------------------------------------------------
2422 subroutine random_init(seed)
2424 C Initialize random number generator
2426 implicit real*8 (a-h,o-z)
2427 include 'DIMENSIONS'
2430 logical OKRandom, prng_restart
2432 integer iseed_array(4)
2434 include 'COMMON.IOUNITS'
2435 include 'COMMON.TIME1'
2436 include 'COMMON.THREAD'
2437 include 'COMMON.SBRIDGE'
2438 include 'COMMON.CONTROL'
2439 include 'COMMON.MCM'
2440 include 'COMMON.MAP'
2441 include 'COMMON.HEADER'
2442 include 'COMMON.CSA'
2443 include 'COMMON.CHAIN'
2444 include 'COMMON.MUCA'
2446 include 'COMMON.FFIELD'
2447 include 'COMMON.SETUP'
2448 iseed=-dint(dabs(seed))
2449 if (iseed.eq.0) then
2450 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2451 & 'Random seed undefined. The program will stop.'
2452 write (*,'(/80(1h*)/20x,a/80(1h*))')
2453 & 'Random seed undefined. The program will stop.'
2455 call mpi_finalize(mpi_comm_world,ierr)
2457 stop 'Bad random seed.'
2460 if (fg_rank.eq.0) then
2464 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2465 OKRandom = prng_restart(me,iseed)
2468 tmp=65536.0d0**(4-i)
2469 iseed_array(i) = dint(seed/tmp)
2470 seed=seed-iseed_array(i)*tmp
2473 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2474 & (iseed_array(i),i=1,4)
2475 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2476 & (iseed_array(i),i=1,4)
2477 OKRandom = prng_restart(me,iseed_array)
2480 c r1 = prng_next(me)
2481 r1=ran_number(0.0D0,1.0D0)
2483 & write (iout,*) 'ran_num',r1
2484 if (r1.lt.0.0d0) OKRandom=.false.
2486 if (.not.OKRandom) then
2487 write (iout,*) 'PRNG IS NOT WORKING!!!'
2488 print *,'PRNG IS NOT WORKING!!!'
2491 call mpi_abort(mpi_comm_world,error_msg,ierr)
2494 write (iout,*) 'too many processors for parallel prng'
2495 write (*,*) 'too many processors for parallel prng'
2503 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)