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 write(iout,*) "bordliptop=",bordliptop
241 write(iout,*) "bordlipbot=",bordlipbot
242 write(iout,*) "bufliptop=",bufliptop
243 write(iout,*) "buflipbot=",buflipbot
246 if (me.eq.king .or. .not.out1file )
247 & write (iout,*) "DISTCHAINMAX",distchainmax
249 if(me.eq.king.or..not.out1file)
250 & write (iout,'(2a)') diagmeth(kdiag),
251 & ' routine used to diagonalize matrices.'
254 c--------------------------------------------------------------------------
255 subroutine read_REMDpar
259 implicit real*8 (a-h,o-z)
261 include 'COMMON.IOUNITS'
262 include 'COMMON.TIME1'
265 include 'COMMON.LANGEVIN'
267 include 'COMMON.LANGEVIN.lang0'
269 include 'COMMON.INTERACT'
270 include 'COMMON.NAMES'
272 include 'COMMON.REMD'
273 include 'COMMON.CONTROL'
274 include 'COMMON.SETUP'
276 character*320 controlcard
277 character*3200 controlcard1
278 integer iremd_m_total
280 if(me.eq.king.or..not.out1file)
281 & write (iout,*) "REMD setup"
283 call card_concat(controlcard)
284 call readi(controlcard,"NREP",nrep,3)
285 call readi(controlcard,"NSTEX",nstex,1000)
286 call reada(controlcard,"RETMIN",retmin,10.0d0)
287 call reada(controlcard,"RETMAX",retmax,1000.0d0)
288 mremdsync=(index(controlcard,'SYNC').gt.0)
289 call readi(controlcard,"NSYN",i_sync_step,100)
290 restart1file=(index(controlcard,'REST1FILE').gt.0)
291 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
292 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
293 if(max_cache_traj_use.gt.max_cache_traj)
294 & max_cache_traj_use=max_cache_traj
295 if(me.eq.king.or..not.out1file) then
296 cd if (traj1file) then
297 crc caching is in testing - NTWX is not ignored
298 cd write (iout,*) "NTWX value is ignored"
299 cd write (iout,*) " trajectory is stored to one file by master"
300 cd write (iout,*) " before exchange at NSTEX intervals"
302 write (iout,*) "NREP= ",nrep
303 write (iout,*) "NSTEX= ",nstex
304 write (iout,*) "SYNC= ",mremdsync
305 write (iout,*) "NSYN= ",i_sync_step
306 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
309 if (index(controlcard,'TLIST').gt.0) then
311 call card_concat(controlcard1)
312 read(controlcard1,*) (remd_t(i),i=1,nrep)
313 if(me.eq.king.or..not.out1file)
314 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
317 if (index(controlcard,'MLIST').gt.0) then
319 call card_concat(controlcard1)
320 read(controlcard1,*) (remd_m(i),i=1,nrep)
321 if(me.eq.king.or..not.out1file) then
322 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
325 iremd_m_total=iremd_m_total+remd_m(i)
327 write (iout,*) 'Total number of replicas ',iremd_m_total
330 if(me.eq.king.or..not.out1file)
331 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
334 c--------------------------------------------------------------------------
335 subroutine read_MDpar
339 implicit real*8 (a-h,o-z)
341 include 'COMMON.IOUNITS'
342 include 'COMMON.TIME1'
345 include 'COMMON.LANGEVIN'
347 include 'COMMON.LANGEVIN.lang0'
349 include 'COMMON.INTERACT'
350 include 'COMMON.NAMES'
352 include 'COMMON.SETUP'
353 include 'COMMON.CONTROL'
354 include 'COMMON.SPLITELE'
356 character*320 controlcard
358 call card_concat(controlcard)
359 call readi(controlcard,"NSTEP",n_timestep,1000000)
360 call readi(controlcard,"NTWE",ntwe,100)
361 call readi(controlcard,"NTWX",ntwx,1000)
362 call reada(controlcard,"DT",d_time,1.0d-1)
363 call reada(controlcard,"DVMAX",dvmax,2.0d1)
364 call reada(controlcard,"DAMAX",damax,1.0d1)
365 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
366 call readi(controlcard,"LANG",lang,0)
367 RESPA = index(controlcard,"RESPA") .gt. 0
368 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
369 ntime_split0=ntime_split
370 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
371 ntime_split0=ntime_split
372 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
373 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
374 rest = index(controlcard,"REST").gt.0
375 tbf = index(controlcard,"TBF").gt.0
376 usampl = index(controlcard,"USAMPL").gt.0
378 mdpdb = index(controlcard,"MDPDB").gt.0
379 call reada(controlcard,"T_BATH",t_bath,300.0d0)
380 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
381 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
382 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
383 if (count_reset_moment.eq.0) count_reset_moment=1000000000
384 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
385 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
386 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
387 if (count_reset_vel.eq.0) count_reset_vel=1000000000
388 large = index(controlcard,"LARGE").gt.0
389 print_compon = index(controlcard,"PRINT_COMPON").gt.0
390 rattle = index(controlcard,"RATTLE").gt.0
391 c if performing umbrella sampling, fragments constrained are read from the fragment file
397 if(me.eq.king.or..not.out1file) then
399 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
401 write (iout,'(a)') "The units are:"
402 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
403 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
404 & " acceleration: angstrom/(48.9 fs)**2"
405 write (iout,'(a)') "energy: kcal/mol, temperature: K"
407 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
408 write (iout,'(a60,f10.5,a)')
409 & "Initial time step of numerical integration:",d_time,
411 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
413 write (iout,'(2a,i4,a)')
414 & "A-MTS algorithm used; initial time step for fast-varying",
415 & " short-range forces split into",ntime_split," steps."
416 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
417 & r_cut," lambda",rlamb
419 write (iout,'(2a,f10.5)')
420 & "Maximum acceleration threshold to reduce the time step",
421 & "/increase split number:",damax
422 write (iout,'(2a,f10.5)')
423 & "Maximum predicted energy drift to reduce the timestep",
424 & "/increase split number:",edriftmax
425 write (iout,'(a60,f10.5)')
426 & "Maximum velocity threshold to reduce velocities:",dvmax
427 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
428 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
429 if (rattle) write (iout,'(a60)')
430 & "Rattle algorithm used to constrain the virtual bonds"
434 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
435 call reada(controlcard,"RWAT",rwat,1.4d0)
436 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
437 surfarea=index(controlcard,"SURFAREA").gt.0
438 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
439 if(me.eq.king.or..not.out1file)then
440 write (iout,'(/a,$)') "Langevin dynamics calculation"
443 & " with direct integration of Langevin equations"
444 else if (lang.eq.2) then
445 write (iout,'(a/)') " with TINKER stochasic MD integrator"
446 else if (lang.eq.3) then
447 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
448 else if (lang.eq.4) then
449 write (iout,'(a/)') " in overdamped mode"
451 write (iout,'(//a,i5)')
452 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
455 write (iout,'(a60,f10.5)') "Temperature:",t_bath
456 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
457 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
458 write (iout,'(a60,f10.5)')
459 & "Scaling factor of the friction forces:",scal_fric
460 if (surfarea) write (iout,'(2a,i10,a)')
461 & "Friction coefficients will be scaled by solvent-accessible",
462 & " surface area every",reset_fricmat," steps."
464 c Calculate friction coefficients and bounds of stochastic forces
465 eta=6*pi*cPoise*etawat
466 if(me.eq.king.or..not.out1file)
467 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
469 gamp=scal_fric*(pstok+rwat)*eta
470 stdfp=dsqrt(2*Rb*t_bath/d_time)
472 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
473 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
475 if(me.eq.king.or..not.out1file)then
476 write (iout,'(/2a/)')
477 & "Radii of site types and friction coefficients and std's of",
478 & " stochastic forces of fully exposed sites"
479 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
481 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
482 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
486 if(me.eq.king.or..not.out1file)then
487 write (iout,'(a)') "Berendsen bath calculation"
488 write (iout,'(a60,f10.5)') "Temperature:",t_bath
489 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
491 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
492 & count_reset_moment," steps"
494 & write (iout,'(a,i10,a)')
495 & "Velocities will be reset at random every",count_reset_vel,
499 if(me.eq.king.or..not.out1file)
500 & write (iout,'(a31)') "Microcanonical mode calculation"
502 if(me.eq.king.or..not.out1file)then
503 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
505 write(iout,*) "MD running with constraints."
506 write(iout,*) "Equilibration time ", eq_time, " mtus."
507 write(iout,*) "Constraining ", nfrag," fragments."
508 write(iout,*) "Length of each fragment, weight and q0:"
510 write (iout,*) "Set of restraints #",iset
512 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
513 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
515 write(iout,*) "constraints between ", npair, "fragments."
516 write(iout,*) "constraint pairs, weights and q0:"
518 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
519 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
521 write(iout,*) "angle constraints within ", nfrag_back,
522 & "backbone fragments."
523 write(iout,*) "fragment, weights:"
525 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
526 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
527 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
530 iset=mod(kolor,nset)+1
533 if(me.eq.king.or..not.out1file)
534 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
537 c------------------------------------------------------------------------------
540 C Read molecular data.
542 implicit real*8 (a-h,o-z)
548 include 'COMMON.IOUNITS'
551 include 'COMMON.INTERACT'
552 include 'COMMON.LOCAL'
553 include 'COMMON.NAMES'
554 include 'COMMON.CHAIN'
555 include 'COMMON.FFIELD'
556 include 'COMMON.SBRIDGE'
557 include 'COMMON.HEADER'
558 include 'COMMON.CONTROL'
559 include 'COMMON.DBASE'
560 include 'COMMON.THREAD'
561 include 'COMMON.CONTACTS'
562 include 'COMMON.TORCNSTR'
563 include 'COMMON.TIME1'
564 include 'COMMON.BOUNDS'
566 include 'COMMON.SETUP'
567 character*4 sequence(maxres)
569 double precision x(maxvar)
570 character*256 pdbfile
571 character*320 weightcard
572 character*80 weightcard_t,ucase
573 dimension itype_pdb(maxres)
574 common /pizda/ itype_pdb
575 logical seq_comp,fail
576 double precision energia(0:n_ene)
582 C Read weights of the subsequent energy terms.
583 call card_concat(weightcard)
584 call reada(weightcard,'WLONG',wlong,1.0D0)
585 call reada(weightcard,'WSC',wsc,wlong)
586 call reada(weightcard,'WSCP',wscp,wlong)
587 call reada(weightcard,'WELEC',welec,1.0D0)
588 call reada(weightcard,'WVDWPP',wvdwpp,welec)
589 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
590 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
591 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
592 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
593 call reada(weightcard,'WTURN3',wturn3,1.0D0)
594 call reada(weightcard,'WTURN4',wturn4,1.0D0)
595 call reada(weightcard,'WTURN6',wturn6,1.0D0)
596 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
597 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
598 call reada(weightcard,'WBOND',wbond,1.0D0)
599 call reada(weightcard,'WTOR',wtor,1.0D0)
600 call reada(weightcard,'WTORD',wtor_d,1.0D0)
601 call reada(weightcard,'WANG',wang,1.0D0)
602 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
603 call reada(weightcard,'SCAL14',scal14,0.4D0)
604 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
605 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
606 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
607 call reada(weightcard,'TEMP0',temp0,300.0d0)
608 call reada(weightcard,'WLT',wliptran,0.0D0)
609 if (index(weightcard,'SOFT').gt.0) ipot=6
610 C 12/1/95 Added weight for the multi-body term WCORR
611 call reada(weightcard,'WCORRH',wcorr,1.0D0)
612 if (wcorr4.gt.0.0d0) wcorr=wcorr4
632 if(me.eq.king.or..not.out1file)
633 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
634 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
636 10 format (/'Energy-term weights (unscaled):'//
637 & 'WSCC= ',f10.6,' (SC-SC)'/
638 & 'WSCP= ',f10.6,' (SC-p)'/
639 & 'WELEC= ',f10.6,' (p-p electr)'/
640 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
641 & 'WBOND= ',f10.6,' (stretching)'/
642 & 'WANG= ',f10.6,' (bending)'/
643 & 'WSCLOC= ',f10.6,' (SC local)'/
644 & 'WTOR= ',f10.6,' (torsional)'/
645 & 'WTORD= ',f10.6,' (double torsional)'/
646 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
647 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
648 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
649 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
650 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
651 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
652 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
653 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
654 & 'WTURN6= ',f10.6,' (turns, 6th order)')
655 if(me.eq.king.or..not.out1file)then
656 if (wcorr4.gt.0.0d0) then
657 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
658 & 'between contact pairs of peptide groups'
659 write (iout,'(2(a,f5.3/))')
660 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
661 & 'Range of quenching the correlation terms:',2*delt_corr
662 else if (wcorr.gt.0.0d0) then
663 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
664 & 'between contact pairs of peptide groups'
666 write (iout,'(a,f8.3)')
667 & 'Scaling factor of 1,4 SC-p interactions:',scal14
668 write (iout,'(a,f8.3)')
669 & 'General scaling factor of SC-p interactions:',scalscp
671 r0_corr=cutoff_corr-delt_corr
673 aad(i,1)=scalscp*aad(i,1)
674 aad(i,2)=scalscp*aad(i,2)
675 bad(i,1)=scalscp*bad(i,1)
676 bad(i,2)=scalscp*bad(i,2)
678 call rescale_weights(t_bath)
679 if(me.eq.king.or..not.out1file)
680 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
681 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
683 22 format (/'Energy-term weights (scaled):'//
684 & 'WSCC= ',f10.6,' (SC-SC)'/
685 & 'WSCP= ',f10.6,' (SC-p)'/
686 & 'WELEC= ',f10.6,' (p-p electr)'/
687 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
688 & 'WBOND= ',f10.6,' (stretching)'/
689 & 'WANG= ',f10.6,' (bending)'/
690 & 'WSCLOC= ',f10.6,' (SC local)'/
691 & 'WTOR= ',f10.6,' (torsional)'/
692 & 'WTORD= ',f10.6,' (double torsional)'/
693 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
694 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
695 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
696 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
697 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
698 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
699 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
700 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
701 & 'WTURN6= ',f10.6,' (turns, 6th order)')
702 if(me.eq.king.or..not.out1file)
703 & write (iout,*) "Reference temperature for weights calculation:",
705 call reada(weightcard,"D0CM",d0cm,3.78d0)
706 call reada(weightcard,"AKCM",akcm,15.1d0)
707 call reada(weightcard,"AKTH",akth,11.0d0)
708 call reada(weightcard,"AKCT",akct,12.0d0)
709 call reada(weightcard,"V1SS",v1ss,-1.08d0)
710 call reada(weightcard,"V2SS",v2ss,7.61d0)
711 call reada(weightcard,"V3SS",v3ss,13.7d0)
712 call reada(weightcard,"EBR",ebr,-5.50D0)
713 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
715 dyn_ss_mask(i)=.false.
719 dyn_ssbond_ij(i,j)=1.0d300
722 call reada(weightcard,"HT",Ht,0.0D0)
724 ss_depth=ebr/wsc-0.25*eps(1,1)
725 Ht=Ht/wsc-0.25*eps(1,1)
726 akcm=akcm*wstrain/wsc
727 akth=akth*wstrain/wsc
728 akct=akct*wstrain/wsc
729 v1ss=v1ss*wstrain/wsc
730 v2ss=v2ss*wstrain/wsc
731 v3ss=v3ss*wstrain/wsc
733 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
736 if(me.eq.king.or..not.out1file) then
737 write (iout,*) "Parameters of the SS-bond potential:"
738 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
740 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
741 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
742 write (iout,*)" HT",Ht
743 print *,'indpdb=',indpdb,' pdbref=',pdbref
745 if (indpdb.gt.0 .or. pdbref) then
746 read(inp,'(a)') pdbfile
747 if(me.eq.king.or..not.out1file)
748 & write (iout,'(2a)') 'PDB data will be read from file ',
749 & pdbfile(:ilen(pdbfile))
750 open(ipdbin,file=pdbfile,status='old',err=33)
752 33 write (iout,'(a)') 'Error opening PDB file.'
755 c print *,'Begin reading pdb data'
757 c print *,'Finished reading pdb data'
758 if(me.eq.king.or..not.out1file)
759 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
760 & ' nstart_sup=',nstart_sup
762 itype_pdb(i)=itype(i)
766 nct=nstart_sup+nsup-1
767 call contact(.false.,ncont_ref,icont_ref,co)
770 if(me.eq.king.or..not.out1file)
771 & write(iout,*)'Adding sidechains'
775 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
778 do while (fail.and.nsi.le.maxsi)
779 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
782 if(fail) write(iout,*)'Adding sidechain failed for res ',
783 & i,' after ',nsi,' trials'
788 if (indpdb.eq.0) then
789 C Read sequence if not taken from the pdb file.
791 c print *,'nres=',nres
792 if (iscode.gt.0) then
793 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
795 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
797 C Convert sequence to numeric code
799 itype(i)=rescode(i,sequence(i),iscode)
801 C Assign initial virtual bond lengths
807 vbld(i+nres)=dsc(iabs(itype(i)))
808 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
809 c write (iout,*) "i",i," itype",itype(i),
810 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
814 c print '(20i4)',(itype(i),i=1,nres)
817 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
819 if (itype(i).eq.ntyp1) then
823 else if (iabs(itype(i+1)).ne.20) then
825 else if (iabs(itype(i)).ne.20) then
832 if(me.eq.king.or..not.out1file)then
833 write (iout,*) "ITEL"
835 write (iout,*) i,itype(i),itel(i)
837 print *,'Call Read_Bridge.'
840 C 8/13/98 Set limits to generating the dihedral angles
845 read (inp,*) ndih_constr
846 if (ndih_constr.gt.0) then
848 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
849 if(me.eq.king.or..not.out1file)then
851 & 'There are',ndih_constr,' constraints on phi angles.'
853 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
857 phi0(i)=deg2rad*phi0(i)
858 drange(i)=deg2rad*drange(i)
860 if(me.eq.king.or..not.out1file)
861 & write (iout,*) 'FTORS',ftors
864 phibound(1,ii) = phi0(i)-drange(i)
865 phibound(2,ii) = phi0(i)+drange(i)
872 write (iout,'(a)') 'Boundaries in phi angle sampling:'
874 write (iout,'(a3,i5,2f10.1)')
875 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
881 cd print *,'NNT=',NNT,' NCT=',NCT
882 if (itype(1).eq.ntyp1) nnt=2
883 if (itype(nres).eq.ntyp1) nct=nct-1
885 if(me.eq.king.or..not.out1file)
886 & write (iout,'(a,i3)') 'nsup=',nsup
888 if (nsup.le.(nct-nnt+1)) then
889 do i=0,nct-nnt+1-nsup
890 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
896 & 'Error - sequences to be superposed do not match.'
899 do i=0,nsup-(nct-nnt+1)
900 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
902 nstart_sup=nstart_sup+i
908 & 'Error - sequences to be superposed do not match.'
911 if (nsup.eq.0) nsup=nct-nnt
912 if (nstart_sup.eq.0) nstart_sup=nnt
913 if (nstart_seq.eq.0) nstart_seq=nnt
914 if(me.eq.king.or..not.out1file)
915 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
916 & ' nstart_seq=',nstart_seq
918 c--- Zscore rms -------
919 if (nz_start.eq.0) nz_start=nnt
920 if (nz_end.eq.0 .and. nsup.gt.0) then
922 else if (nz_end.eq.0) then
925 if(me.eq.king.or..not.out1file)then
926 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
927 write (iout,*) 'IZ_SC=',iz_sc
929 c----------------------
932 if (.not.pdbref) then
933 call read_angles(inp,*38)
935 38 write (iout,'(a)') 'Error reading reference structure.'
937 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
938 stop 'Error reading reference structure'
942 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
952 call contact(.true.,ncont_ref,icont_ref,co)
954 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
956 if (constr_dist.gt.0) call read_dist_constr
957 write (iout,*) "After read_dist_constr nhpb",nhpb
959 if(me.eq.king.or..not.out1file)
960 & write (iout,*) 'Contact order:',co
962 if(me.eq.king.or..not.out1file)
963 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
966 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
968 if(me.eq.king.or..not.out1file)
969 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
970 & icont_ref(1,i),' ',
971 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
975 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
976 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
977 & modecalc.ne.10) then
978 C If input structure hasn't been supplied from the PDB file read or generate
980 if (iranconf.eq.0 .and. .not. extconf) then
981 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
982 & write (iout,'(a)') 'Initial geometry will be read in.'
984 read(inp,'(8f10.5)',end=36,err=36)
985 & ((c(l,k),l=1,3),k=1,nres),
986 & ((c(l,k+nres),l=1,3),k=nnt,nct)
987 write (iout,*) "Exit READ_CART"
988 write (iout,'(8f10.5)')
989 & ((c(l,k),l=1,3),k=1,nres),
990 & ((c(l,k+nres),l=1,3),k=nnt,nct)
991 call int_from_cart1(.true.)
992 write (iout,*) "Finish INT_TO_CART"
995 dc(j,i)=c(j,i+1)-c(j,i)
996 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1000 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1002 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1003 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1009 call read_angles(inp,*36)
1012 36 write (iout,'(a)') 'Error reading angle file.'
1014 call mpi_finalize( MPI_COMM_WORLD,IERR )
1016 stop 'Error reading angle file.'
1018 else if (extconf) then
1019 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1020 & write (iout,'(a)') 'Extended chain initial geometry.'
1022 theta(i)=90d0*deg2rad
1025 phi(i)=180d0*deg2rad
1028 alph(i)=110d0*deg2rad
1031 omeg(i)=-120d0*deg2rad
1032 if (itype(i).le.0) omeg(i)=-omeg(i)
1035 if(me.eq.king.or..not.out1file)
1036 & write (iout,'(a)') 'Random-generated initial geometry.'
1040 if (me.eq.king .or. fg_rank.eq.0 .and. (
1041 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1045 call gen_rand_conf(itmp,*30)
1047 30 write (iout,*) 'Failed to generate random conformation',
1048 & ', itrial=',itrial
1049 write (*,*) 'Processor:',me,
1050 & ' Failed to generate random conformation',
1060 write (iout,'(a,i3,a)') 'Processor:',me,
1061 & ' error in generating random conformation.'
1062 write (*,'(a,i3,a)') 'Processor:',me,
1063 & ' error in generating random conformation.'
1066 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1071 & ' error in generating random conformation.'
1076 elseif (modecalc.eq.4) then
1077 read (inp,'(a)') intinname
1078 open (intin,file=intinname,status='old',err=333)
1079 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1080 & write (iout,'(a)') 'intinname',intinname
1081 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1083 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1085 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1087 stop 'Error opening angle file.'
1091 C Generate distance constraints, if the PDB structure is to be regularized.
1092 if (nthread.gt.0) then
1093 call read_threadbase
1096 if (me.eq.king .or. .not. out1file)
1098 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1099 write (iout,'(/a,i3,a)')
1100 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1101 write (iout,'(20i4)') (iss(i),i=1,ns)
1103 write(iout,*)"Running with dynamic disulfide-bond formation"
1105 write (iout,'(/a/)') 'Pre-formed links are:'
1111 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1112 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1118 if (ns.gt.0.and.dyn_ss) then
1122 forcon(i-nss)=forcon(i)
1129 dyn_ss_mask(iss(i))=.true.
1132 if (i2ndstr.gt.0) call secstrp2dihc
1133 c call geom_to_var(nvar,x)
1134 c call etotal(energia(0))
1135 c call enerprint(energia(0))
1136 c call briefout(0,etot)
1138 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1139 cd write (iout,'(a)') 'Variable list:'
1140 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1142 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1143 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1144 & 'Processor',myrank,': end reading molecular data.'
1148 c--------------------------------------------------------------------------
1149 logical function seq_comp(itypea,itypeb,length)
1151 integer length,itypea(length),itypeb(length)
1154 if (itypea(i).ne.itypeb(i)) then
1162 c-----------------------------------------------------------------------------
1163 subroutine read_bridge
1164 C Read information about disulfide bridges.
1165 implicit real*8 (a-h,o-z)
1166 include 'DIMENSIONS'
1170 include 'COMMON.IOUNITS'
1171 include 'COMMON.GEO'
1172 include 'COMMON.VAR'
1173 include 'COMMON.INTERACT'
1174 include 'COMMON.LOCAL'
1175 include 'COMMON.NAMES'
1176 include 'COMMON.CHAIN'
1177 include 'COMMON.FFIELD'
1178 include 'COMMON.SBRIDGE'
1179 include 'COMMON.HEADER'
1180 include 'COMMON.CONTROL'
1181 include 'COMMON.DBASE'
1182 include 'COMMON.THREAD'
1183 include 'COMMON.TIME1'
1184 include 'COMMON.SETUP'
1185 C Read bridging residues.
1186 read (inp,*) ns,(iss(i),i=1,ns)
1188 if(me.eq.king.or..not.out1file)
1189 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1190 C Check whether the specified bridging residues are cystines.
1192 if (itype(iss(i)).ne.1) then
1193 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1194 & 'Do you REALLY think that the residue ',
1195 & restyp(itype(iss(i))),i,
1196 & ' can form a disulfide bridge?!!!'
1197 write (*,'(2a,i3,a)')
1198 & 'Do you REALLY think that the residue ',
1199 & restyp(itype(iss(i))),i,
1200 & ' can form a disulfide bridge?!!!'
1202 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1207 C Read preformed bridges.
1209 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1211 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1214 C Check if the residues involved in bridges are in the specified list of
1215 C bridging residues.
1218 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1219 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1220 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1221 & ' contains residues present in other pairs.'
1222 write (*,'(a,i3,a)') 'Disulfide pair',i,
1223 & ' contains residues present in other pairs.'
1225 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1231 if (ihpb(i).eq.iss(j)) goto 10
1233 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1236 if (jhpb(i).eq.iss(j)) goto 20
1238 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1244 ihpb(i)=ihpb(i)+nres
1245 jhpb(i)=jhpb(i)+nres
1251 c----------------------------------------------------------------------------
1252 subroutine read_x(kanal,*)
1253 implicit real*8 (a-h,o-z)
1254 include 'DIMENSIONS'
1255 include 'COMMON.GEO'
1256 include 'COMMON.VAR'
1257 include 'COMMON.CHAIN'
1258 include 'COMMON.IOUNITS'
1259 include 'COMMON.CONTROL'
1260 include 'COMMON.LOCAL'
1261 include 'COMMON.INTERACT'
1262 c Read coordinates from input
1264 read(kanal,'(8f10.5)',end=10,err=10)
1265 & ((c(l,k),l=1,3),k=1,nres),
1266 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1269 c(j,2*nres)=c(j,nres)
1271 call int_from_cart1(.false.)
1274 dc(j,i)=c(j,i+1)-c(j,i)
1275 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1279 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1281 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1282 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1290 c----------------------------------------------------------------------------
1291 subroutine read_threadbase
1292 implicit real*8 (a-h,o-z)
1293 include 'DIMENSIONS'
1294 include 'COMMON.IOUNITS'
1295 include 'COMMON.GEO'
1296 include 'COMMON.VAR'
1297 include 'COMMON.INTERACT'
1298 include 'COMMON.LOCAL'
1299 include 'COMMON.NAMES'
1300 include 'COMMON.CHAIN'
1301 include 'COMMON.FFIELD'
1302 include 'COMMON.SBRIDGE'
1303 include 'COMMON.HEADER'
1304 include 'COMMON.CONTROL'
1305 include 'COMMON.DBASE'
1306 include 'COMMON.THREAD'
1307 include 'COMMON.TIME1'
1308 C Read pattern database for threading.
1309 read (icbase,*) nseq
1311 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1312 & nres_base(2,i),nres_base(3,i)
1313 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1315 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1316 c & nres_base(2,i),nres_base(3,i)
1317 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1321 if (weidis.eq.0.0D0) weidis=0.1D0
1330 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1331 write (iout,'(a,i5)') 'nexcl: ',nexcl
1332 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1335 c------------------------------------------------------------------------------
1336 subroutine setup_var
1337 implicit real*8 (a-h,o-z)
1338 include 'DIMENSIONS'
1339 include 'COMMON.IOUNITS'
1340 include 'COMMON.GEO'
1341 include 'COMMON.VAR'
1342 include 'COMMON.INTERACT'
1343 include 'COMMON.LOCAL'
1344 include 'COMMON.NAMES'
1345 include 'COMMON.CHAIN'
1346 include 'COMMON.FFIELD'
1347 include 'COMMON.SBRIDGE'
1348 include 'COMMON.HEADER'
1349 include 'COMMON.CONTROL'
1350 include 'COMMON.DBASE'
1351 include 'COMMON.THREAD'
1352 include 'COMMON.TIME1'
1353 C Set up variable list.
1359 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1361 ialph(i,1)=nvar+nside
1365 if (indphi.gt.0) then
1367 else if (indback.gt.0) then
1372 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1375 c----------------------------------------------------------------------------
1376 subroutine gen_dist_constr
1377 C Generate CA distance constraints.
1378 implicit real*8 (a-h,o-z)
1379 include 'DIMENSIONS'
1380 include 'COMMON.IOUNITS'
1381 include 'COMMON.GEO'
1382 include 'COMMON.VAR'
1383 include 'COMMON.INTERACT'
1384 include 'COMMON.LOCAL'
1385 include 'COMMON.NAMES'
1386 include 'COMMON.CHAIN'
1387 include 'COMMON.FFIELD'
1388 include 'COMMON.SBRIDGE'
1389 include 'COMMON.HEADER'
1390 include 'COMMON.CONTROL'
1391 include 'COMMON.DBASE'
1392 include 'COMMON.THREAD'
1393 include 'COMMON.TIME1'
1394 dimension itype_pdb(maxres)
1395 common /pizda/ itype_pdb
1397 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1398 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1399 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1401 do i=nstart_sup,nstart_sup+nsup-1
1402 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1403 cd & ' seq_pdb', restyp(itype_pdb(i))
1404 do j=i+2,nstart_sup+nsup-1
1406 ihpb(nhpb)=i+nstart_seq-nstart_sup
1407 jhpb(nhpb)=j+nstart_seq-nstart_sup
1409 dhpb(nhpb)=dist(i,j)
1412 cd write (iout,'(a)') 'Distance constraints:'
1417 cd if (ii.gt.nres) then
1422 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1423 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1424 cd & dhpb(i),forcon(i)
1428 c----------------------------------------------------------------------------
1430 implicit real*8 (a-h,o-z)
1431 include 'DIMENSIONS'
1432 include 'COMMON.MAP'
1433 include 'COMMON.IOUNITS'
1434 character*3 angid(4) /'THE','PHI','ALP','OME'/
1435 character*80 mapcard,ucase
1437 read (inp,'(a)') mapcard
1438 mapcard=ucase(mapcard)
1439 if (index(mapcard,'PHI').gt.0) then
1441 else if (index(mapcard,'THE').gt.0) then
1443 else if (index(mapcard,'ALP').gt.0) then
1445 else if (index(mapcard,'OME').gt.0) then
1448 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1449 stop 'Error - illegal variable spec in MAP card.'
1451 call readi (mapcard,'RES1',res1(imap),0)
1452 call readi (mapcard,'RES2',res2(imap),0)
1453 if (res1(imap).eq.0) then
1454 res1(imap)=res2(imap)
1455 else if (res2(imap).eq.0) then
1456 res2(imap)=res1(imap)
1458 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1460 & 'Error - illegal definition of variable group in MAP.'
1461 stop 'Error - illegal definition of variable group in MAP.'
1463 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1464 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1465 call readi(mapcard,'NSTEP',nstep(imap),0)
1466 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1468 & 'Illegal boundary and/or step size specification in MAP.'
1469 stop 'Illegal boundary and/or step size specification in MAP.'
1474 c----------------------------------------------------------------------------
1476 implicit real*8 (a-h,o-z)
1477 include 'DIMENSIONS'
1478 include 'COMMON.IOUNITS'
1479 include 'COMMON.GEO'
1480 include 'COMMON.CSA'
1481 include 'COMMON.BANK'
1482 include 'COMMON.CONTROL'
1484 character*620 mcmcard
1485 call card_concat(mcmcard)
1487 call readi(mcmcard,'NCONF',nconf,50)
1488 call readi(mcmcard,'NADD',nadd,0)
1489 call readi(mcmcard,'JSTART',jstart,1)
1490 call readi(mcmcard,'JEND',jend,1)
1491 call readi(mcmcard,'NSTMAX',nstmax,500000)
1492 call readi(mcmcard,'N0',n0,1)
1493 call readi(mcmcard,'N1',n1,6)
1494 call readi(mcmcard,'N2',n2,4)
1495 call readi(mcmcard,'N3',n3,0)
1496 call readi(mcmcard,'N4',n4,0)
1497 call readi(mcmcard,'N5',n5,0)
1498 call readi(mcmcard,'N6',n6,10)
1499 call readi(mcmcard,'N7',n7,0)
1500 call readi(mcmcard,'N8',n8,0)
1501 call readi(mcmcard,'N9',n9,0)
1502 call readi(mcmcard,'N14',n14,0)
1503 call readi(mcmcard,'N15',n15,0)
1504 call readi(mcmcard,'N16',n16,0)
1505 call readi(mcmcard,'N17',n17,0)
1506 call readi(mcmcard,'N18',n18,0)
1508 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1510 call readi(mcmcard,'NDIFF',ndiff,2)
1511 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1512 call readi(mcmcard,'IS1',is1,1)
1513 call readi(mcmcard,'IS2',is2,8)
1514 call readi(mcmcard,'NRAN0',nran0,4)
1515 call readi(mcmcard,'NRAN1',nran1,2)
1516 call readi(mcmcard,'IRR',irr,1)
1517 call readi(mcmcard,'NSEED',nseed,20)
1518 call readi(mcmcard,'NTOTAL',ntotal,10000)
1519 call reada(mcmcard,'CUT1',cut1,2.0d0)
1520 call reada(mcmcard,'CUT2',cut2,5.0d0)
1521 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1522 call readi(mcmcard,'ICMAX',icmax,3)
1523 call readi(mcmcard,'IRESTART',irestart,0)
1524 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1527 call reada(mcmcard,'DELE',dele,20.0d0)
1528 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1529 call readi(mcmcard,'IREF',iref,0)
1530 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1531 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1532 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1533 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1534 write (iout,*) "NCONF_IN",nconf_in
1537 c----------------------------------------------------------------------------
1538 cfmc subroutine mcmfread
1539 cfmc implicit real*8 (a-h,o-z)
1540 cfmc include 'DIMENSIONS'
1541 cfmc include 'COMMON.MCMF'
1542 cfmc include 'COMMON.IOUNITS'
1543 cfmc include 'COMMON.GEO'
1544 cfmc character*80 ucase
1545 cfmc character*620 mcmcard
1546 cfmc call card_concat(mcmcard)
1548 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1549 cfmc write(iout,*)'MAXRANT=',maxrant
1550 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1551 cfmc write(iout,*)'MAXFAM=',maxfam
1552 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1553 cfmc write(iout,*)'NNET1=',nnet1
1554 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1555 cfmc write(iout,*)'NNET2=',nnet2
1556 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1557 cfmc write(iout,*)'NNET3=',nnet3
1558 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1559 cfmc write(iout,*)'ILASTT=',ilastt
1560 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1561 cfmc write(iout,*)'MAXSTR=',maxstr
1562 cfmc maxstr_f=maxstr/maxfam
1563 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1564 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1565 cfmc write(iout,*)'NMCMF=',nmcmf
1566 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1567 cfmc write(iout,*)'IFOCUS=',ifocus
1568 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1569 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1570 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1571 cfmc write(iout,*)'INTPRT=',intprt
1572 cfmc call readi(mcmcard,'IPRT',iprt,100)
1573 cfmc write(iout,*)'IPRT=',iprt
1574 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1575 cfmc write(iout,*)'IMAXTR=',imaxtr
1576 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1577 cfmc write(iout,*)'MAXEVEN=',maxeven
1578 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1579 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1580 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1581 cfmc write(iout,*)'INIMIN=',inimin
1582 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1583 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1584 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1585 cfmc write(iout,*)'NTHREAD=',nthread
1586 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1587 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1588 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1589 cfmc write(iout,*)'MAXPERT=',maxpert
1590 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1591 cfmc write(iout,*)'IRMSD=',irmsd
1592 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1593 cfmc write(iout,*)'DENEMIN=',denemin
1594 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1595 cfmc write(iout,*)'RCUT1S=',rcut1s
1596 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1597 cfmc write(iout,*)'RCUT1E=',rcut1e
1598 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1599 cfmc write(iout,*)'RCUT2S=',rcut2s
1600 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1601 cfmc write(iout,*)'RCUT2E=',rcut2e
1602 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1603 cfmc write(iout,*)'DPERT1=',d_pert1
1604 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1605 cfmc write(iout,*)'DPERT1A=',d_pert1a
1606 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1607 cfmc write(iout,*)'DPERT2=',d_pert2
1608 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1609 cfmc write(iout,*)'DPERT2A=',d_pert2a
1610 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1611 cfmc write(iout,*)'DPERT2B=',d_pert2b
1612 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1613 cfmc write(iout,*)'DPERT2C=',d_pert2c
1614 cfmc d_pert1=deg2rad*d_pert1
1615 cfmc d_pert1a=deg2rad*d_pert1a
1616 cfmc d_pert2=deg2rad*d_pert2
1617 cfmc d_pert2a=deg2rad*d_pert2a
1618 cfmc d_pert2b=deg2rad*d_pert2b
1619 cfmc d_pert2c=deg2rad*d_pert2c
1620 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1621 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1622 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1623 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1624 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1625 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1626 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1627 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1628 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1629 cfmc write(iout,*)'RCUTINI=',rcutini
1630 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1631 cfmc write(iout,*)'GRAT=',grat
1632 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1633 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1637 c----------------------------------------------------------------------------
1639 implicit real*8 (a-h,o-z)
1640 include 'DIMENSIONS'
1641 include 'COMMON.MCM'
1642 include 'COMMON.MCE'
1643 include 'COMMON.IOUNITS'
1645 character*320 mcmcard
1646 call card_concat(mcmcard)
1647 call readi(mcmcard,'MAXACC',maxacc,100)
1648 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1649 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1650 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1651 call readi(mcmcard,'MAXREPM',maxrepm,200)
1652 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1653 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1654 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1655 call reada(mcmcard,'E_UP',e_up,5.0D0)
1656 call reada(mcmcard,'DELTE',delte,0.1D0)
1657 call readi(mcmcard,'NSWEEP',nsweep,5)
1658 call readi(mcmcard,'NSTEPH',nsteph,0)
1659 call readi(mcmcard,'NSTEPC',nstepc,0)
1660 call reada(mcmcard,'TMIN',tmin,298.0D0)
1661 call reada(mcmcard,'TMAX',tmax,298.0D0)
1662 call readi(mcmcard,'NWINDOW',nwindow,0)
1663 call readi(mcmcard,'PRINT_MC',print_mc,0)
1664 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1665 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1666 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1667 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1668 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1669 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1670 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1671 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1672 if (nwindow.gt.0) then
1673 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1675 winlen(i)=winend(i)-winstart(i)+1
1678 if (tmax.lt.tmin) tmax=tmin
1679 if (tmax.eq.tmin) then
1683 if (nstepc.gt.0 .and. nsteph.gt.0) then
1684 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1685 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1687 C Probabilities of different move types
1688 sumpro_type(0)=0.0D0
1689 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1690 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1691 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1692 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1693 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1694 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1695 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1697 print *,'i',i,' sumprotype',sumpro_type(i)
1698 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1699 print *,'i',i,' sumprotype',sumpro_type(i)
1703 c----------------------------------------------------------------------------
1704 subroutine read_minim
1705 implicit real*8 (a-h,o-z)
1706 include 'DIMENSIONS'
1707 include 'COMMON.MINIM'
1708 include 'COMMON.IOUNITS'
1710 character*320 minimcard
1711 call card_concat(minimcard)
1712 call readi(minimcard,'MAXMIN',maxmin,2000)
1713 call readi(minimcard,'MAXFUN',maxfun,5000)
1714 call readi(minimcard,'MINMIN',minmin,maxmin)
1715 call readi(minimcard,'MINFUN',minfun,maxmin)
1716 call reada(minimcard,'TOLF',tolf,1.0D-2)
1717 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1718 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1719 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1720 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1721 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1722 & 'Options in energy minimization:'
1723 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1724 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1725 & 'MinMin:',MinMin,' MinFun:',MinFun,
1726 & ' TolF:',TolF,' RTolF:',RTolF
1729 c----------------------------------------------------------------------------
1730 subroutine read_angles(kanal,*)
1731 implicit real*8 (a-h,o-z)
1732 include 'DIMENSIONS'
1733 include 'COMMON.GEO'
1734 include 'COMMON.VAR'
1735 include 'COMMON.CHAIN'
1736 include 'COMMON.IOUNITS'
1737 include 'COMMON.CONTROL'
1738 c Read angles from input
1740 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1741 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1742 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1743 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1746 c 9/7/01 avoid 180 deg valence angle
1747 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1749 theta(i)=deg2rad*theta(i)
1750 phi(i)=deg2rad*phi(i)
1751 alph(i)=deg2rad*alph(i)
1752 omeg(i)=deg2rad*omeg(i)
1757 c----------------------------------------------------------------------------
1758 subroutine reada(rekord,lancuch,wartosc,default)
1760 character*(*) rekord,lancuch
1761 double precision wartosc,default
1764 iread=index(rekord,lancuch)
1765 if (iread.eq.0) then
1769 iread=iread+ilen(lancuch)+1
1770 read (rekord(iread:),*,err=10,end=10) wartosc
1775 c----------------------------------------------------------------------------
1776 subroutine readi(rekord,lancuch,wartosc,default)
1778 character*(*) rekord,lancuch
1779 integer wartosc,default
1782 iread=index(rekord,lancuch)
1783 if (iread.eq.0) then
1787 iread=iread+ilen(lancuch)+1
1788 read (rekord(iread:),*,err=10,end=10) wartosc
1793 c----------------------------------------------------------------------------
1794 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1797 integer tablica(dim),default
1798 character*(*) rekord,lancuch
1805 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1806 if (iread.eq.0) return
1807 iread=iread+ilen(lancuch)+1
1808 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1811 c----------------------------------------------------------------------------
1812 subroutine multreada(rekord,lancuch,tablica,dim,default)
1815 double precision tablica(dim),default
1816 character*(*) rekord,lancuch
1823 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1824 if (iread.eq.0) return
1825 iread=iread+ilen(lancuch)+1
1826 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1829 c----------------------------------------------------------------------------
1830 subroutine openunits
1831 implicit real*8 (a-h,o-z)
1832 include 'DIMENSIONS'
1835 character*16 form,nodename
1838 include 'COMMON.SETUP'
1839 include 'COMMON.IOUNITS'
1841 include 'COMMON.CONTROL'
1842 integer lenpre,lenpot,ilen,lentmp
1844 character*3 out1file_text,ucase
1847 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1848 call getenv_loc("PREFIX",prefix)
1850 call getenv_loc("POT",pot)
1851 call getenv_loc("DIRTMP",tmpdir)
1852 call getenv_loc("CURDIR",curdir)
1853 call getenv_loc("OUT1FILE",out1file_text)
1854 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1855 out1file_text=ucase(out1file_text)
1856 if (out1file_text(1:1).eq."Y") then
1859 out1file=fg_rank.gt.0
1864 if (lentmp.gt.0) then
1865 write (*,'(80(1h!))')
1866 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
1867 write (*,'(80(1h!))')
1868 write (*,*)"All output files will be on node /tmp directory."
1870 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1871 if (me.eq.king) then
1872 write (*,*) "The master node is ",nodename
1873 else if (fg_rank.eq.0) then
1874 write (*,*) "I am the CG slave node ",nodename
1876 write (*,*) "I am the FG slave node ",nodename
1879 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1880 lenpre = lentmp+lenpre+1
1882 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1883 C Get the names and open the input files
1884 #if defined(WINIFL) || defined(WINPGI)
1885 open(1,file=pref_orig(:ilen(pref_orig))//
1886 & '.inp',status='old',readonly,shared)
1887 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1888 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1889 C Get parameter filenames and open the parameter files.
1890 call getenv_loc('BONDPAR',bondname)
1891 open (ibond,file=bondname,status='old',readonly,shared)
1892 call getenv_loc('THETPAR',thetname)
1893 open (ithep,file=thetname,status='old',readonly,shared)
1894 call getenv_loc('ROTPAR',rotname)
1895 open (irotam,file=rotname,status='old',readonly,shared)
1896 call getenv_loc('TORPAR',torname)
1897 open (itorp,file=torname,status='old',readonly,shared)
1898 call getenv_loc('TORDPAR',tordname)
1899 open (itordp,file=tordname,status='old',readonly,shared)
1900 call getenv_loc('FOURIER',fouriername)
1901 open (ifourier,file=fouriername,status='old',readonly,shared)
1902 call getenv_loc('ELEPAR',elename)
1903 open (ielep,file=elename,status='old',readonly,shared)
1904 call getenv_loc('SIDEPAR',sidename)
1905 open (isidep,file=sidename,status='old',readonly,shared)
1906 #elif (defined CRAY) || (defined AIX)
1907 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1909 c print *,"Processor",myrank," opened file 1"
1910 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1911 c print *,"Processor",myrank," opened file 9"
1912 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1913 C Get parameter filenames and open the parameter files.
1914 call getenv_loc('BONDPAR',bondname)
1915 open (ibond,file=bondname,status='old',action='read')
1916 c print *,"Processor",myrank," opened file IBOND"
1917 call getenv_loc('THETPAR',thetname)
1918 open (ithep,file=thetname,status='old',action='read')
1919 c print *,"Processor",myrank," opened file ITHEP"
1920 call getenv_loc('ROTPAR',rotname)
1921 open (irotam,file=rotname,status='old',action='read')
1922 c print *,"Processor",myrank," opened file IROTAM"
1923 call getenv_loc('TORPAR',torname)
1924 open (itorp,file=torname,status='old',action='read')
1925 c print *,"Processor",myrank," opened file ITORP"
1926 call getenv_loc('TORDPAR',tordname)
1927 open (itordp,file=tordname,status='old',action='read')
1928 c print *,"Processor",myrank," opened file ITORDP"
1929 call getenv_loc('SCCORPAR',sccorname)
1930 open (isccor,file=sccorname,status='old',action='read')
1931 c print *,"Processor",myrank," opened file ISCCOR"
1932 call getenv_loc('FOURIER',fouriername)
1933 open (ifourier,file=fouriername,status='old',action='read')
1934 c print *,"Processor",myrank," opened file IFOURIER"
1935 call getenv_loc('ELEPAR',elename)
1936 open (ielep,file=elename,status='old',action='read')
1937 c print *,"Processor",myrank," opened file IELEP"
1938 call getenv_loc('SIDEPAR',sidename)
1939 open (isidep,file=sidename,status='old',action='read')
1940 c print *,"Processor",myrank," opened file ISIDEP"
1941 c print *,"Processor",myrank," opened parameter files"
1943 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1944 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1945 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1946 C Get parameter filenames and open the parameter files.
1947 call getenv_loc('BONDPAR',bondname)
1948 open (ibond,file=bondname,status='old')
1949 call getenv_loc('THETPAR',thetname)
1950 open (ithep,file=thetname,status='old')
1951 call getenv_loc('ROTPAR',rotname)
1952 open (irotam,file=rotname,status='old')
1953 call getenv_loc('TORPAR',torname)
1954 open (itorp,file=torname,status='old')
1955 call getenv_loc('TORDPAR',tordname)
1956 open (itordp,file=tordname,status='old')
1957 call getenv_loc('SCCORPAR',sccorname)
1958 open (isccor,file=sccorname,status='old')
1959 call getenv_loc('FOURIER',fouriername)
1960 open (ifourier,file=fouriername,status='old')
1961 call getenv_loc('ELEPAR',elename)
1962 open (ielep,file=elename,status='old')
1963 call getenv_loc('SIDEPAR',sidename)
1964 open (isidep,file=sidename,status='old')
1966 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1968 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1969 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1970 C Get parameter filenames and open the parameter files.
1971 call getenv_loc('BONDPAR',bondname)
1972 open (ibond,file=bondname,status='old',readonly)
1973 call getenv_loc('THETPAR',thetname)
1974 open (ithep,file=thetname,status='old',readonly)
1975 call getenv_loc('ROTPAR',rotname)
1976 open (irotam,file=rotname,status='old',readonly)
1977 call getenv_loc('TORPAR',torname)
1978 open (itorp,file=torname,status='old',readonly)
1979 call getenv_loc('TORDPAR',tordname)
1980 open (itordp,file=tordname,status='old',readonly)
1981 call getenv_loc('SCCORPAR',sccorname)
1982 open (isccor,file=sccorname,status='old',readonly)
1984 call getenv_loc('THETPARPDB',thetname_pdb)
1985 print *,"thetname_pdb ",thetname_pdb
1986 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1987 print *,ithep_pdb," opened"
1989 call getenv_loc('FOURIER',fouriername)
1990 open (ifourier,file=fouriername,status='old',readonly)
1991 call getenv_loc('ELEPAR',elename)
1992 open (ielep,file=elename,status='old',readonly)
1993 call getenv_loc('SIDEPAR',sidename)
1994 open (isidep,file=sidename,status='old',readonly)
1995 call getenv_loc('LIPTRANPAR',liptranname)
1996 open (iliptranpar,file=liptranname,status='old',action='read')
1998 call getenv_loc('ROTPARPDB',rotname_pdb)
1999 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2004 C 8/9/01 In the newest version SCp interaction constants are read from a file
2005 C Use -DOLDSCP to use hard-coded constants instead.
2007 call getenv_loc('SCPPAR',scpname)
2008 #if defined(WINIFL) || defined(WINPGI)
2009 open (iscpp,file=scpname,status='old',readonly,shared)
2010 #elif (defined CRAY) || (defined AIX)
2011 open (iscpp,file=scpname,status='old',action='read')
2013 open (iscpp,file=scpname,status='old')
2015 open (iscpp,file=scpname,status='old',readonly)
2018 call getenv_loc('PATTERN',patname)
2019 #if defined(WINIFL) || defined(WINPGI)
2020 open (icbase,file=patname,status='old',readonly,shared)
2021 #elif (defined CRAY) || (defined AIX)
2022 open (icbase,file=patname,status='old',action='read')
2024 open (icbase,file=patname,status='old')
2026 open (icbase,file=patname,status='old',readonly)
2029 C Open output file only for CG processes
2030 c print *,"Processor",myrank," fg_rank",fg_rank
2031 if (fg_rank.eq.0) then
2033 if (nodes.eq.1) then
2036 npos = dlog10(dfloat(nodes-1))+1
2038 if (npos.lt.3) npos=3
2039 write (liczba,'(i1)') npos
2040 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2042 write (liczba,form) me
2043 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2044 & liczba(:ilen(liczba))
2045 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2047 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2049 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2050 & liczba(:ilen(liczba))//'.mol2'
2051 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2052 & liczba(:ilen(liczba))//'.stat'
2054 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2055 & //liczba(:ilen(liczba))//'.stat')
2056 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2059 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2060 & liczba(:ilen(liczba))//'.const'
2065 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2066 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2067 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2068 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2069 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2071 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2073 rest2name=prefix(:ilen(prefix))//'.rst'
2075 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2078 #if defined(AIX) || defined(PGI)
2079 if (me.eq.king .or. .not. out1file)
2080 & open(iout,file=outname,status='unknown')
2082 if (fg_rank.gt.0) then
2083 write (liczba,'(i3.3)') myrank/nfgtasks
2084 write (ll,'(bz,i3.3)') fg_rank
2085 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2090 open(igeom,file=intname,status='unknown',position='append')
2091 open(ipdb,file=pdbname,status='unknown')
2092 open(imol2,file=mol2name,status='unknown')
2093 open(istat,file=statname,status='unknown',position='append')
2095 c1out open(iout,file=outname,status='unknown')
2098 if (me.eq.king .or. .not.out1file)
2099 & open(iout,file=outname,status='unknown')
2101 if (fg_rank.gt.0) then
2102 write (liczba,'(i3.3)') myrank/nfgtasks
2103 write (ll,'(bz,i3.3)') fg_rank
2104 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2109 open(igeom,file=intname,status='unknown',access='append')
2110 open(ipdb,file=pdbname,status='unknown')
2111 open(imol2,file=mol2name,status='unknown')
2112 open(istat,file=statname,status='unknown',access='append')
2114 c1out open(iout,file=outname,status='unknown')
2117 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2118 csa_seed=prefix(:lenpre)//'.CSA.seed'
2119 csa_history=prefix(:lenpre)//'.CSA.history'
2120 csa_bank=prefix(:lenpre)//'.CSA.bank'
2121 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2122 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2123 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2124 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2125 csa_int=prefix(:lenpre)//'.int'
2126 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2127 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2128 csa_in=prefix(:lenpre)//'.CSA.in'
2129 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2132 write (iout,'(80(1h-))')
2133 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2134 write (iout,'(80(1h-))')
2135 write (iout,*) "Input file : ",
2136 & pref_orig(:ilen(pref_orig))//'.inp'
2137 write (iout,*) "Output file : ",
2138 & outname(:ilen(outname))
2140 write (iout,*) "Sidechain potential file : ",
2141 & sidename(:ilen(sidename))
2143 write (iout,*) "SCp potential file : ",
2144 & scpname(:ilen(scpname))
2146 write (iout,*) "Electrostatic potential file : ",
2147 & elename(:ilen(elename))
2148 write (iout,*) "Cumulant coefficient file : ",
2149 & fouriername(:ilen(fouriername))
2150 write (iout,*) "Torsional parameter file : ",
2151 & torname(:ilen(torname))
2152 write (iout,*) "Double torsional parameter file : ",
2153 & tordname(:ilen(tordname))
2154 write (iout,*) "SCCOR parameter file : ",
2155 & sccorname(:ilen(sccorname))
2156 write (iout,*) "Bond & inertia constant file : ",
2157 & bondname(:ilen(bondname))
2158 write (iout,*) "Bending parameter file : ",
2159 & thetname(:ilen(thetname))
2160 write (iout,*) "Rotamer parameter file : ",
2161 & rotname(:ilen(rotname))
2162 write (iout,*) "Thetpdb parameter file : ",
2163 & thetname_pdb(:ilen(thetname_pdb))
2164 write (iout,*) "Threading database : ",
2165 & patname(:ilen(patname))
2167 &write (iout,*)" DIRTMP : ",
2169 write (iout,'(80(1h-))')
2173 c----------------------------------------------------------------------------
2174 subroutine card_concat(card)
2175 implicit real*8 (a-h,o-z)
2176 include 'DIMENSIONS'
2177 include 'COMMON.IOUNITS'
2179 character*80 karta,ucase
2181 read (inp,'(a)') karta
2184 do while (karta(80:80).eq.'&')
2185 card=card(:ilen(card)+1)//karta(:79)
2186 read (inp,'(a)') karta
2189 card=card(:ilen(card)+1)//karta
2192 c----------------------------------------------------------------------------------
2194 implicit real*8 (a-h,o-z)
2195 include 'DIMENSIONS'
2196 include 'COMMON.CHAIN'
2197 include 'COMMON.IOUNITS'
2199 open(irest2,file=rest2name,status='unknown')
2200 read(irest2,*) totT,EK,potE,totE,t_bath
2202 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2205 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2208 read (irest2,*) iset
2213 c---------------------------------------------------------------------------------
2214 subroutine read_fragments
2215 implicit real*8 (a-h,o-z)
2216 include 'DIMENSIONS'
2220 include 'COMMON.SETUP'
2221 include 'COMMON.CHAIN'
2222 include 'COMMON.IOUNITS'
2224 include 'COMMON.CONTROL'
2225 read(inp,*) nset,nfrag,npair,nfrag_back
2226 if(me.eq.king.or..not.out1file)
2227 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2228 & " nfrag_back",nfrag_back
2230 read(inp,*) mset(iset)
2232 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2234 if(me.eq.king.or..not.out1file)
2235 & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2236 & ifrag(2,i,iset), qinfrag(i,iset)
2239 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2241 if(me.eq.king.or..not.out1file)
2242 & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2243 & ipair(2,i,iset), qinpair(i,iset)
2246 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2247 & wfrag_back(3,i,iset),
2248 & ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2249 if(me.eq.king.or..not.out1file)
2250 & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2251 & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2256 c-------------------------------------------------------------------------------
2257 subroutine read_dist_constr
2258 implicit real*8 (a-h,o-z)
2259 include 'DIMENSIONS'
2263 include 'COMMON.SETUP'
2264 include 'COMMON.CONTROL'
2265 include 'COMMON.CHAIN'
2266 include 'COMMON.IOUNITS'
2267 include 'COMMON.SBRIDGE'
2268 integer ifrag_(2,100),ipair_(2,100)
2269 double precision wfrag_(100),wpair_(100)
2270 character*500 controlcard
2271 c write (iout,*) "Calling read_dist_constr"
2272 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2274 call card_concat(controlcard)
2275 call readi(controlcard,"NFRAG",nfrag_,0)
2276 call readi(controlcard,"NPAIR",npair_,0)
2277 call readi(controlcard,"NDIST",ndist_,0)
2278 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2279 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2280 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2281 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2282 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2283 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2284 c write (iout,*) "IFRAG"
2286 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2288 c write (iout,*) "IPAIR"
2290 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2294 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2295 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2296 & ifrag_(2,i)=nstart_sup+nsup-1
2297 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2299 if (wfrag_(i).gt.0.0d0) then
2300 do j=ifrag_(1,i),ifrag_(2,i)-1
2301 do k=j+1,ifrag_(2,i)
2302 c write (iout,*) "j",j," k",k
2304 if (constr_dist.eq.1) then
2309 forcon(nhpb)=wfrag_(i)
2310 else if (constr_dist.eq.2) then
2311 if (ddjk.le.dist_cut) then
2316 forcon(nhpb)=wfrag_(i)
2323 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2326 if (.not.out1file .or. me.eq.king)
2327 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2328 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2330 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2331 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2338 if (wpair_(i).gt.0.0d0) then
2346 do j=ifrag_(1,ii),ifrag_(2,ii)
2347 do k=ifrag_(1,jj),ifrag_(2,jj)
2351 forcon(nhpb)=wpair_(i)
2352 dhpb(nhpb)=dist(j,k)
2354 if (.not.out1file .or. me.eq.king)
2355 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2356 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2358 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2359 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2366 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2367 if (forcon(nhpb+1).gt.0.0d0) then
2369 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2371 if (.not.out1file .or. me.eq.king)
2372 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2373 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2375 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2376 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2383 c-------------------------------------------------------------------------------
2385 subroutine flush(iu)
2390 subroutine flush(iu)
2395 c------------------------------------------------------------------------------
2396 subroutine copy_to_tmp(source)
2397 include "DIMENSIONS"
2398 include "COMMON.IOUNITS"
2399 character*(*) source
2400 character* 256 tmpfile
2404 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2405 inquire(file=tmpfile,exist=ex)
2407 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2408 & " to temporary directory..."
2409 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2410 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2414 c------------------------------------------------------------------------------
2415 subroutine move_from_tmp(source)
2416 include "DIMENSIONS"
2417 include "COMMON.IOUNITS"
2418 character*(*) source
2421 write (*,*) "Moving ",source(:ilen(source)),
2422 & " from temporary directory to working directory"
2423 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2424 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2427 c------------------------------------------------------------------------------
2428 subroutine random_init(seed)
2430 C Initialize random number generator
2432 implicit real*8 (a-h,o-z)
2433 include 'DIMENSIONS'
2436 logical OKRandom, prng_restart
2438 integer iseed_array(4)
2440 include 'COMMON.IOUNITS'
2441 include 'COMMON.TIME1'
2442 include 'COMMON.THREAD'
2443 include 'COMMON.SBRIDGE'
2444 include 'COMMON.CONTROL'
2445 include 'COMMON.MCM'
2446 include 'COMMON.MAP'
2447 include 'COMMON.HEADER'
2448 include 'COMMON.CSA'
2449 include 'COMMON.CHAIN'
2450 include 'COMMON.MUCA'
2452 include 'COMMON.FFIELD'
2453 include 'COMMON.SETUP'
2454 iseed=-dint(dabs(seed))
2455 if (iseed.eq.0) then
2456 write (iout,'(/80(1h*)/20x,a/80(1h*))')
2457 & 'Random seed undefined. The program will stop.'
2458 write (*,'(/80(1h*)/20x,a/80(1h*))')
2459 & 'Random seed undefined. The program will stop.'
2461 call mpi_finalize(mpi_comm_world,ierr)
2463 stop 'Bad random seed.'
2466 if (fg_rank.eq.0) then
2470 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2471 OKRandom = prng_restart(me,iseed)
2474 tmp=65536.0d0**(4-i)
2475 iseed_array(i) = dint(seed/tmp)
2476 seed=seed-iseed_array(i)*tmp
2479 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2480 & (iseed_array(i),i=1,4)
2481 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2482 & (iseed_array(i),i=1,4)
2483 OKRandom = prng_restart(me,iseed_array)
2486 c r1 = prng_next(me)
2487 r1=ran_number(0.0D0,1.0D0)
2489 & write (iout,*) 'ran_num',r1
2490 if (r1.lt.0.0d0) OKRandom=.false.
2492 if (.not.OKRandom) then
2493 write (iout,*) 'PRNG IS NOT WORKING!!!'
2494 print *,'PRNG IS NOT WORKING!!!'
2497 call mpi_abort(mpi_comm_world,error_msg,ierr)
2500 write (iout,*) 'too many processors for parallel prng'
2501 write (*,*) 'too many processors for parallel prng'
2509 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)