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 write (iout,*) "constr_dist",constr_dist
102 call readi(controlcard,'NSAXS',nsaxs,0)
103 call readi(controlcard,'SAXS_MODE',saxs_mode,0)
104 call reada(controlcard,'SCAL_RAD',scal_rad,1.0d0)
105 call reada(controlcard,'SAXS_CUTOFF',saxs_cutoff,1.0d0)
106 write (iout,*) "Number of SAXS restraints",NSAXS," SAXS_MODE",
107 & SAXS_MODE," SCAL_RAD",scal_rad,"SAXS_CUTOFF",saxs_cutoff
108 call readi(controlcard,'CONSTR_HOMOL',constr_homology,0)
109 call readi(controlcard,'SYM',symetr,1)
110 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
111 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
112 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
113 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
114 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
115 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
116 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
117 call reada(controlcard,'DRMS',drms,0.1D0)
118 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
119 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
120 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
121 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
122 write (iout,'(a,f10.1)')'DRMS = ',drms
123 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
124 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
126 call readi(controlcard,'NZ_START',nz_start,0)
127 call readi(controlcard,'NZ_END',nz_end,0)
128 call readi(controlcard,'IZ_SC',iz_sc,0)
130 safety = 60.0d0*safety
133 call reada(controlcard,"T_BATH",t_bath,300.0d0)
134 minim=(index(controlcard,'MINIMIZE').gt.0)
135 dccart=(index(controlcard,'CART').gt.0)
136 overlapsc=(index(controlcard,'OVERLAP').gt.0)
137 overlapsc=.not.overlapsc
138 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
139 searchsc=.not.searchsc
140 sideadd=(index(controlcard,'SIDEADD').gt.0)
141 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
142 outpdb=(index(controlcard,'PDBOUT').gt.0)
143 outmol2=(index(controlcard,'MOL2OUT').gt.0)
144 pdbref=(index(controlcard,'PDBREF').gt.0)
145 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
146 indpdb=index(controlcard,'PDBSTART')
147 extconf=(index(controlcard,'EXTCONF').gt.0)
148 AFMlog=(index(controlcard,'AFM'))
149 selfguide=(index(controlcard,'SELFGUIDE'))
150 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
151 call readi(controlcard,'IPRINT',iprint,0)
152 call readi(controlcard,'MAXGEN',maxgen,10000)
153 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
154 call readi(controlcard,"KDIAG",kdiag,0)
155 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
156 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
157 & write (iout,*) "RESCALE_MODE",rescale_mode
158 split_ene=index(controlcard,'SPLIT_ENE').gt.0
159 if (index(controlcard,'REGULAR').gt.0.0D0) then
160 call reada(controlcard,'WEIDIS',weidis,0.1D0)
164 if (index(controlcard,'CHECKGRAD').gt.0) then
166 if (index(controlcard,'CART').gt.0) then
168 elseif (index(controlcard,'CARINT').gt.0) then
173 elseif (index(controlcard,'THREAD').gt.0) then
175 call readi(controlcard,'THREAD',nthread,0)
176 if (nthread.gt.0) then
177 call reada(controlcard,'WEIDIS',weidis,0.1D0)
180 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
181 stop 'Error termination in Read_Control.'
183 else if (index(controlcard,'MCMA').gt.0) then
185 else if (index(controlcard,'MCEE').gt.0) then
187 else if (index(controlcard,'MULTCONF').gt.0) then
189 else if (index(controlcard,'MAP').gt.0) then
191 call readi(controlcard,'MAP',nmap,0)
192 else if (index(controlcard,'CSA').gt.0) then
194 crc else if (index(controlcard,'ZSCORE').gt.0) then
196 crc ZSCORE is rm from UNRES, modecalc=9 is available
199 cfcm else if (index(controlcard,'MCMF').gt.0) then
201 else if (index(controlcard,'SOFTREG').gt.0) then
203 else if (index(controlcard,'CHECK_BOND').gt.0) then
205 else if (index(controlcard,'TEST').gt.0) then
207 else if (index(controlcard,'MD').gt.0) then
209 else if (index(controlcard,'RE ').gt.0) then
213 lmuca=index(controlcard,'MUCA').gt.0
214 call readi(controlcard,'MUCADYN',mucadyn,0)
215 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
216 if (lmuca .and. (me.eq.king .or. .not.out1file ))
218 write (iout,*) 'MUCADYN=',mucadyn
219 write (iout,*) 'MUCASMOOTH=',muca_smooth
222 iscode=index(controlcard,'ONE_LETTER')
223 indphi=index(controlcard,'PHI')
224 indback=index(controlcard,'BACK')
225 iranconf=index(controlcard,'RAND_CONF')
226 i2ndstr=index(controlcard,'USE_SEC_PRED')
227 gradout=index(controlcard,'GRADOUT').gt.0
228 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
229 C DISTCHAINMAX become obsolete for periodic boundry condition
230 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
231 C Reading the dimensions of box in x,y,z coordinates
232 call reada(controlcard,'BOXX',boxxsize,100.0d0)
233 call reada(controlcard,'BOXY',boxysize,100.0d0)
234 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
235 c Cutoff range for interactions
236 call reada(controlcard,"R_CUT",r_cut,15.0d0)
237 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
238 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
239 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
240 if (lipthick.gt.0.0d0) then
241 bordliptop=(boxzsize+lipthick)/2.0
242 bordlipbot=bordliptop-lipthick
244 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
245 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
246 buflipbot=bordlipbot+lipbufthick
247 bufliptop=bordliptop-lipbufthick
248 if ((lipbufthick*2.0d0).gt.lipthick)
249 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
251 c write(iout,*) "bordliptop=",bordliptop
252 c write(iout,*) "bordlipbot=",bordlipbot
253 c write(iout,*) "bufliptop=",bufliptop
254 c write(iout,*) "buflipbot=",buflipbot
257 if (me.eq.king .or. .not.out1file )
258 & write (iout,*) "DISTCHAINMAX",distchainmax
260 if(me.eq.king.or..not.out1file)
261 & write (iout,'(2a)') diagmeth(kdiag),
262 & ' routine used to diagonalize matrices.'
265 c--------------------------------------------------------------------------
266 subroutine read_REMDpar
270 implicit real*8 (a-h,o-z)
272 include 'COMMON.IOUNITS'
273 include 'COMMON.TIME1'
276 include 'COMMON.LANGEVIN'
278 include 'COMMON.LANGEVIN.lang0'
280 include 'COMMON.INTERACT'
281 include 'COMMON.NAMES'
283 include 'COMMON.REMD'
284 include 'COMMON.CONTROL'
285 include 'COMMON.SETUP'
287 character*320 controlcard
288 character*3200 controlcard1
289 integer iremd_m_total
291 if(me.eq.king.or..not.out1file)
292 & write (iout,*) "REMD setup"
294 call card_concat(controlcard)
295 call readi(controlcard,"NREP",nrep,3)
296 call readi(controlcard,"NSTEX",nstex,1000)
297 call reada(controlcard,"RETMIN",retmin,10.0d0)
298 call reada(controlcard,"RETMAX",retmax,1000.0d0)
299 mremdsync=(index(controlcard,'SYNC').gt.0)
300 call readi(controlcard,"NSYN",i_sync_step,100)
301 restart1file=(index(controlcard,'REST1FILE').gt.0)
302 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
303 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
304 if(max_cache_traj_use.gt.max_cache_traj)
305 & max_cache_traj_use=max_cache_traj
306 if(me.eq.king.or..not.out1file) then
307 cd if (traj1file) then
308 crc caching is in testing - NTWX is not ignored
309 cd write (iout,*) "NTWX value is ignored"
310 cd write (iout,*) " trajectory is stored to one file by master"
311 cd write (iout,*) " before exchange at NSTEX intervals"
313 write (iout,*) "NREP= ",nrep
314 write (iout,*) "NSTEX= ",nstex
315 write (iout,*) "SYNC= ",mremdsync
316 write (iout,*) "NSYN= ",i_sync_step
317 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
320 if (index(controlcard,'TLIST').gt.0) then
322 call card_concat(controlcard1)
323 read(controlcard1,*) (remd_t(i),i=1,nrep)
324 if(me.eq.king.or..not.out1file)
325 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
328 if (index(controlcard,'MLIST').gt.0) then
330 call card_concat(controlcard1)
331 read(controlcard1,*) (remd_m(i),i=1,nrep)
332 if(me.eq.king.or..not.out1file) then
333 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
336 iremd_m_total=iremd_m_total+remd_m(i)
338 write (iout,*) 'Total number of replicas ',iremd_m_total
341 if(me.eq.king.or..not.out1file)
342 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
345 c--------------------------------------------------------------------------
346 subroutine read_MDpar
350 implicit real*8 (a-h,o-z)
352 include 'COMMON.IOUNITS'
353 include 'COMMON.TIME1'
356 include 'COMMON.LANGEVIN'
358 include 'COMMON.LANGEVIN.lang0'
360 include 'COMMON.INTERACT'
361 include 'COMMON.NAMES'
363 include 'COMMON.SETUP'
364 include 'COMMON.CONTROL'
365 include 'COMMON.SPLITELE'
367 character*320 controlcard
369 call card_concat(controlcard)
370 call readi(controlcard,"NSTEP",n_timestep,1000000)
371 call readi(controlcard,"NTWE",ntwe,100)
372 call readi(controlcard,"NTWX",ntwx,1000)
373 call reada(controlcard,"DT",d_time,1.0d-1)
374 call reada(controlcard,"DVMAX",dvmax,2.0d1)
375 call reada(controlcard,"DAMAX",damax,1.0d1)
376 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
377 call readi(controlcard,"LANG",lang,0)
378 RESPA = index(controlcard,"RESPA") .gt. 0
379 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
380 ntime_split0=ntime_split
381 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
382 ntime_split0=ntime_split
383 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
384 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
385 rest = index(controlcard,"REST").gt.0
386 tbf = index(controlcard,"TBF").gt.0
387 usampl = index(controlcard,"USAMPL").gt.0
389 mdpdb = index(controlcard,"MDPDB").gt.0
390 call reada(controlcard,"T_BATH",t_bath,300.0d0)
391 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
392 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
393 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
394 if (count_reset_moment.eq.0) count_reset_moment=1000000000
395 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
396 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
397 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
398 if (count_reset_vel.eq.0) count_reset_vel=1000000000
399 large = index(controlcard,"LARGE").gt.0
400 print_compon = index(controlcard,"PRINT_COMPON").gt.0
401 rattle = index(controlcard,"RATTLE").gt.0
402 preminim = index(controlcard,"PREMINIM").gt.0
404 dccart=(index(controlcard,'CART').gt.0)
407 c if performing umbrella sampling, fragments constrained are read from the fragment file
413 if(me.eq.king.or..not.out1file) then
415 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
417 write (iout,'(a)') "The units are:"
418 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
419 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
420 & " acceleration: angstrom/(48.9 fs)**2"
421 write (iout,'(a)') "energy: kcal/mol, temperature: K"
423 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
424 write (iout,'(a60,f10.5,a)')
425 & "Initial time step of numerical integration:",d_time,
427 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
429 write (iout,'(2a,i4,a)')
430 & "A-MTS algorithm used; initial time step for fast-varying",
431 & " short-range forces split into",ntime_split," steps."
432 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
433 & r_cut," lambda",rlamb
435 write (iout,'(2a,f10.5)')
436 & "Maximum acceleration threshold to reduce the time step",
437 & "/increase split number:",damax
438 write (iout,'(2a,f10.5)')
439 & "Maximum predicted energy drift to reduce the timestep",
440 & "/increase split number:",edriftmax
441 write (iout,'(a60,f10.5)')
442 & "Maximum velocity threshold to reduce velocities:",dvmax
443 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
444 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
445 if (rattle) write (iout,'(a60)')
446 & "Rattle algorithm used to constrain the virtual bonds"
447 if (preminim .or. iranconf.gt.0) then
449 & "Initial structure will be energy-minimized"
454 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
455 call reada(controlcard,"RWAT",rwat,1.4d0)
456 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
457 surfarea=index(controlcard,"SURFAREA").gt.0
458 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
459 if(me.eq.king.or..not.out1file)then
460 write (iout,'(/a,$)') "Langevin dynamics calculation"
463 & " with direct integration of Langevin equations"
464 else if (lang.eq.2) then
465 write (iout,'(a/)') " with TINKER stochasic MD integrator"
466 else if (lang.eq.3) then
467 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
468 else if (lang.eq.4) then
469 write (iout,'(a/)') " in overdamped mode"
471 write (iout,'(//a,i5)')
472 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
475 write (iout,'(a60,f10.5)') "Temperature:",t_bath
476 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
477 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
478 write (iout,'(a60,f10.5)')
479 & "Scaling factor of the friction forces:",scal_fric
480 if (surfarea) write (iout,'(2a,i10,a)')
481 & "Friction coefficients will be scaled by solvent-accessible",
482 & " surface area every",reset_fricmat," steps."
484 c Calculate friction coefficients and bounds of stochastic forces
485 eta=6*pi*cPoise*etawat
486 if(me.eq.king.or..not.out1file)
487 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
489 gamp=scal_fric*(pstok+rwat)*eta
490 stdfp=dsqrt(2*Rb*t_bath/d_time)
492 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
493 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
495 if(me.eq.king.or..not.out1file)then
496 write (iout,'(/2a/)')
497 & "Radii of site types and friction coefficients and std's of",
498 & " stochastic forces of fully exposed sites"
499 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
501 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
502 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
506 if(me.eq.king.or..not.out1file)then
507 write (iout,'(a)') "Berendsen bath calculation"
508 write (iout,'(a60,f10.5)') "Temperature:",t_bath
509 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
511 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
512 & count_reset_moment," steps"
514 & write (iout,'(a,i10,a)')
515 & "Velocities will be reset at random every",count_reset_vel,
519 if(me.eq.king.or..not.out1file)
520 & write (iout,'(a31)') "Microcanonical mode calculation"
522 if(me.eq.king.or..not.out1file)then
523 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
525 write(iout,*) "MD running with constraints."
526 write(iout,*) "Equilibration time ", eq_time, " mtus."
527 write(iout,*) "Constraining ", nfrag," fragments."
528 write(iout,*) "Length of each fragment, weight and q0:"
530 write (iout,*) "Set of restraints #",iset
532 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
533 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
535 write(iout,*) "constraints between ", npair, "fragments."
536 write(iout,*) "constraint pairs, weights and q0:"
538 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
539 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
541 write(iout,*) "angle constraints within ", nfrag_back,
542 & "backbone fragments."
543 write(iout,*) "fragment, weights:"
545 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
546 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
547 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
550 iset=mod(kolor,nset)+1
553 if(me.eq.king.or..not.out1file)
554 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
557 c------------------------------------------------------------------------------
560 C Read molecular data.
562 implicit real*8 (a-h,o-z)
568 include 'COMMON.IOUNITS'
571 include 'COMMON.INTERACT'
572 include 'COMMON.LOCAL'
573 include 'COMMON.NAMES'
574 include 'COMMON.CHAIN'
575 include 'COMMON.FFIELD'
576 include 'COMMON.SBRIDGE'
577 include 'COMMON.HEADER'
578 include 'COMMON.CONTROL'
579 include 'COMMON.DBASE'
580 include 'COMMON.THREAD'
581 include 'COMMON.CONTACTS'
582 include 'COMMON.TORCNSTR'
583 include 'COMMON.TIME1'
584 include 'COMMON.BOUNDS'
586 include 'COMMON.SETUP'
587 character*4 sequence(maxres)
589 double precision x(maxvar)
590 character*256 pdbfile
591 character*320 weightcard
592 character*80 weightcard_t,ucase
593 dimension itype_pdb(maxres)
594 common /pizda/ itype_pdb
595 logical seq_comp,fail
596 double precision energia(0:n_ene)
602 C Read weights of the subsequent energy terms.
603 call card_concat(weightcard)
604 call reada(weightcard,'WLONG',wlong,1.0D0)
605 call reada(weightcard,'WSC',wsc,wlong)
606 call reada(weightcard,'WSCP',wscp,wlong)
607 call reada(weightcard,'WELEC',welec,1.0D0)
608 call reada(weightcard,'WVDWPP',wvdwpp,welec)
609 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
610 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
611 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
612 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
613 call reada(weightcard,'WTURN3',wturn3,1.0D0)
614 call reada(weightcard,'WTURN4',wturn4,1.0D0)
615 call reada(weightcard,'WTURN6',wturn6,1.0D0)
616 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
617 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
618 call reada(weightcard,'WBOND',wbond,1.0D0)
619 call reada(weightcard,'WTOR',wtor,1.0D0)
620 call reada(weightcard,'WTORD',wtor_d,1.0D0)
621 call reada(weightcard,'WANG',wang,1.0D0)
622 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
623 call reada(weightcard,'SCAL14',scal14,0.4D0)
624 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
625 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
626 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
627 call reada(weightcard,'TEMP0',temp0,300.0d0)
628 call reada(weightcard,'WLT',wliptran,0.0D0)
629 call reada(weightcard,'WSAXS',wsaxs,1.0D0)
630 if (index(weightcard,'SOFT').gt.0) ipot=6
631 C 12/1/95 Added weight for the multi-body term WCORR
632 call reada(weightcard,'WCORRH',wcorr,1.0D0)
633 if (wcorr4.gt.0.0d0) wcorr=wcorr4
654 if(me.eq.king.or..not.out1file)
655 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
656 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
658 10 format (/'Energy-term weights (unscaled):'//
659 & 'WSCC= ',f10.6,' (SC-SC)'/
660 & 'WSCP= ',f10.6,' (SC-p)'/
661 & 'WELEC= ',f10.6,' (p-p electr)'/
662 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
663 & 'WBOND= ',f10.6,' (stretching)'/
664 & 'WANG= ',f10.6,' (bending)'/
665 & 'WSCLOC= ',f10.6,' (SC local)'/
666 & 'WTOR= ',f10.6,' (torsional)'/
667 & 'WTORD= ',f10.6,' (double torsional)'/
668 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
669 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
670 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
671 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
672 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
673 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
674 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
675 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
676 & 'WTURN6= ',f10.6,' (turns, 6th order)')
677 if(me.eq.king.or..not.out1file)then
678 if (wcorr4.gt.0.0d0) then
679 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
680 & 'between contact pairs of peptide groups'
681 write (iout,'(2(a,f5.3/))')
682 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
683 & 'Range of quenching the correlation terms:',2*delt_corr
684 else if (wcorr.gt.0.0d0) then
685 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
686 & 'between contact pairs of peptide groups'
688 write (iout,'(a,f8.3)')
689 & 'Scaling factor of 1,4 SC-p interactions:',scal14
690 write (iout,'(a,f8.3)')
691 & 'General scaling factor of SC-p interactions:',scalscp
693 r0_corr=cutoff_corr-delt_corr
695 aad(i,1)=scalscp*aad(i,1)
696 aad(i,2)=scalscp*aad(i,2)
697 bad(i,1)=scalscp*bad(i,1)
698 bad(i,2)=scalscp*bad(i,2)
700 call rescale_weights(t_bath)
701 if(me.eq.king.or..not.out1file)
702 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
703 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
705 22 format (/'Energy-term weights (scaled):'//
706 & 'WSCC= ',f10.6,' (SC-SC)'/
707 & 'WSCP= ',f10.6,' (SC-p)'/
708 & 'WELEC= ',f10.6,' (p-p electr)'/
709 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
710 & 'WBOND= ',f10.6,' (stretching)'/
711 & 'WANG= ',f10.6,' (bending)'/
712 & 'WSCLOC= ',f10.6,' (SC local)'/
713 & 'WTOR= ',f10.6,' (torsional)'/
714 & 'WTORD= ',f10.6,' (double torsional)'/
715 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
716 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
717 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
718 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
719 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
720 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
721 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
722 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
723 & 'WTURN6= ',f10.6,' (turns, 6th order)')
724 if(me.eq.king.or..not.out1file)
725 & write (iout,*) "Reference temperature for weights calculation:",
727 call reada(weightcard,"D0CM",d0cm,3.78d0)
728 call reada(weightcard,"AKCM",akcm,15.1d0)
729 call reada(weightcard,"AKTH",akth,11.0d0)
730 call reada(weightcard,"AKCT",akct,12.0d0)
731 call reada(weightcard,"V1SS",v1ss,-1.08d0)
732 call reada(weightcard,"V2SS",v2ss,7.61d0)
733 call reada(weightcard,"V3SS",v3ss,13.7d0)
734 call reada(weightcard,"EBR",ebr,-5.50D0)
735 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
737 dyn_ss_mask(i)=.false.
741 dyn_ssbond_ij(i,j)=1.0d300
744 call reada(weightcard,"HT",Ht,0.0D0)
746 ss_depth=ebr/wsc-0.25*eps(1,1)
747 Ht=Ht/wsc-0.25*eps(1,1)
748 akcm=akcm*wstrain/wsc
749 akth=akth*wstrain/wsc
750 akct=akct*wstrain/wsc
751 v1ss=v1ss*wstrain/wsc
752 v2ss=v2ss*wstrain/wsc
753 v3ss=v3ss*wstrain/wsc
755 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
758 if(me.eq.king.or..not.out1file) then
759 write (iout,*) "Parameters of the SS-bond potential:"
760 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
762 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
763 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
764 write (iout,*)" HT",Ht
765 print *,'indpdb=',indpdb,' pdbref=',pdbref
767 if (indpdb.gt.0 .or. pdbref) then
768 read(inp,'(a)') pdbfile
769 if(me.eq.king.or..not.out1file)
770 & write (iout,'(2a)') 'PDB data will be read from file ',
771 & pdbfile(:ilen(pdbfile))
772 open(ipdbin,file=pdbfile,status='old',err=33)
774 33 write (iout,'(a)') 'Error opening PDB file.'
777 c print *,'Begin reading pdb data'
786 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
787 & (crefjlee(j,i+nres),j=1,3)
790 c print *,'Finished reading pdb data'
791 if(me.eq.king.or..not.out1file)
792 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
793 & ' nstart_sup=',nstart_sup
795 itype_pdb(i)=itype(i)
799 nct=nstart_sup+nsup-1
800 call contact(.false.,ncont_ref,icont_ref,co)
803 if(me.eq.king.or..not.out1file)
804 & write(iout,*)'Adding sidechains'
808 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
811 do while (fail.and.nsi.le.maxsi)
812 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
816 c Calculalte only the coordinates of the current sidechain; no need to rebuild
818 call locate_side_chain(i)
819 if(fail) write(iout,*)'Adding sidechain failed for res ',
820 & i,' after ',nsi,' trials'
823 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
827 if (indpdb.eq.0) then
828 C Read sequence if not taken from the pdb file.
830 c print *,'nres=',nres
831 if (iscode.gt.0) then
832 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
834 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
836 C Convert sequence to numeric code
838 itype(i)=rescode(i,sequence(i),iscode)
840 C Assign initial virtual bond lengths
846 vbld(i+nres)=dsc(iabs(itype(i)))
847 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
848 c write (iout,*) "i",i," itype",itype(i),
849 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
853 c print '(20i4)',(itype(i),i=1,nres)
856 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
858 if (itype(i).eq.ntyp1) then
862 else if (iabs(itype(i+1)).ne.20) then
864 else if (iabs(itype(i)).ne.20) then
871 if(me.eq.king.or..not.out1file)then
872 write (iout,*) "ITEL"
874 write (iout,*) i,itype(i),itel(i)
876 print *,'Call Read_Bridge.'
879 C 8/13/98 Set limits to generating the dihedral angles
884 read (inp,*) ndih_constr
885 write (iout,*) "ndish_constr",ndih_constr
886 if (ndih_constr.gt.0) then
888 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
889 if(me.eq.king.or..not.out1file)then
891 & 'There are',ndih_constr,' constraints on phi angles.'
893 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
897 phi0(i)=deg2rad*phi0(i)
898 drange(i)=deg2rad*drange(i)
900 if(me.eq.king.or..not.out1file)
901 & write (iout,*) 'FTORS',ftors
904 phibound(1,ii) = phi0(i)-drange(i)
905 phibound(2,ii) = phi0(i)+drange(i)
912 write (iout,'(a)') 'Boundaries in phi angle sampling:'
914 write (iout,'(a3,i5,2f10.1)')
915 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
921 cd print *,'NNT=',NNT,' NCT=',NCT
922 if (itype(1).eq.ntyp1) nnt=2
923 if (itype(nres).eq.ntyp1) nct=nct-1
925 if(me.eq.king.or..not.out1file)
926 & write (iout,'(a,i3)') 'nsup=',nsup
928 if (nsup.le.(nct-nnt+1)) then
929 do i=0,nct-nnt+1-nsup
930 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
936 & 'Error - sequences to be superposed do not match.'
939 do i=0,nsup-(nct-nnt+1)
940 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
942 nstart_sup=nstart_sup+i
948 & 'Error - sequences to be superposed do not match.'
951 if (nsup.eq.0) nsup=nct-nnt
952 if (nstart_sup.eq.0) nstart_sup=nnt
953 if (nstart_seq.eq.0) nstart_seq=nnt
954 if(me.eq.king.or..not.out1file)
955 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
956 & ' nstart_seq=',nstart_seq
958 c--- Zscore rms -------
959 if (nz_start.eq.0) nz_start=nnt
960 if (nz_end.eq.0 .and. nsup.gt.0) then
962 else if (nz_end.eq.0) then
965 if(me.eq.king.or..not.out1file)then
966 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
967 write (iout,*) 'IZ_SC=',iz_sc
969 c----------------------
972 if (.not.pdbref) then
973 call read_angles(inp,*38)
975 38 write (iout,'(a)') 'Error reading reference structure.'
977 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
978 stop 'Error reading reference structure'
980 39 call chainbuild_extconf
982 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
992 call contact(.true.,ncont_ref,icont_ref,co)
994 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
996 if (constr_dist.gt.0) call read_dist_constr
997 c write (iout,*) "After read_dist_constr nhpb",nhpb
998 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
999 if(me.eq.king.or..not.out1file)
1000 & write (iout,*) 'Contact order:',co
1002 if(me.eq.king.or..not.out1file)
1003 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1006 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1008 if(me.eq.king.or..not.out1file)
1009 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1010 & icont_ref(1,i),' ',
1011 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1015 write (iout,*) "calling read_saxs_consrtr",nsaxs
1016 if (nsaxs.gt.0) call read_saxs_constr
1019 if (constr_homology.gt.0) then
1020 call read_constr_homology
1021 if (indpdb.gt.0 .or. pdbref) then
1024 c(j,i)=crefjlee(j,i)
1025 cref(j,i,1)=crefjlee(j,i)
1030 write (iout,*) "Array C"
1032 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1033 & (c(j,i+nres),j=1,3)
1035 write (iout,*) "Array Cref"
1037 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1038 & (cref(j,i+nres,1),j=1,3)
1041 call int_from_cart1(.false.)
1042 call sc_loc_geom(.false.)
1044 thetaref(i)=theta(i)
1049 dc(j,i)=c(j,i+1)-c(j,i)
1050 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1055 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1056 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1064 if (nhpb.gt.0) call hpb_partition
1065 c write (iout,*) "After read_dist_constr nhpb",nhpb
1067 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1068 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1069 & modecalc.ne.10) then
1070 C If input structure hasn't been supplied from the PDB file read or generate
1072 if (iranconf.eq.0 .and. .not. extconf) then
1073 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1074 & write (iout,'(a)') 'Initial geometry will be read in.'
1076 read(inp,'(8f10.5)',end=36,err=36)
1077 & ((c(l,k),l=1,3),k=1,nres),
1078 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1079 write (iout,*) "Exit READ_CART"
1080 write (iout,'(8f10.5)')
1081 & ((c(l,k),l=1,3),k=1,nres),
1082 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1083 call int_from_cart1(.true.)
1084 write (iout,*) "Finish INT_TO_CART"
1087 dc(j,i)=c(j,i+1)-c(j,i)
1088 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1092 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1094 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1095 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1101 write (iout,*) "Calling read_ang"
1102 call read_angles(inp,*36)
1103 write (iout,*) "Calling chainbuild"
1104 call chainbuild_extconf
1107 36 write (iout,'(a)') 'Error reading angle file.'
1109 call mpi_finalize( MPI_COMM_WORLD,IERR )
1111 stop 'Error reading angle file.'
1113 else if (extconf) then
1114 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1115 & write (iout,'(a)') 'Extended chain initial geometry.'
1117 theta(i)=90d0*deg2rad
1120 phi(i)=180d0*deg2rad
1123 alph(i)=110d0*deg2rad
1126 omeg(i)=-120d0*deg2rad
1127 if (itype(i).le.0) omeg(i)=-omeg(i)
1129 c from old chainbuild
1131 C Define the origin and orientation of the coordinate system and locate the
1132 C first three CA's and SC(2).
1136 * Build the alpha-carbon chain.
1139 call locate_next_res(i)
1142 C First and last SC must coincide with the corresponding CA.
1146 dc_norm(j,nres+1)=0.0D0
1147 dc(j,nres+nres)=0.0D0
1148 dc_norm(j,nres+nres)=0.0D0
1150 c(j,nres+nres)=c(j,nres)
1153 C Define the origin and orientation of the coordinate system and locate the
1154 C first three CA's and SC(2).
1158 * Build the alpha-carbon chain.
1161 call locate_next_res(i)
1164 C First and last SC must coincide with the corresponding CA.
1168 dc_norm(j,nres+1)=0.0D0
1169 dc(j,nres+nres)=0.0D0
1170 dc_norm(j,nres+nres)=0.0D0
1172 c(j,nres+nres)=c(j,nres)
1177 if(me.eq.king.or..not.out1file)
1178 & write (iout,'(a)') 'Random-generated initial geometry.'
1182 if (me.eq.king .or. fg_rank.eq.0 .and. (
1183 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1187 call gen_rand_conf(itmp,*30)
1189 30 write (iout,*) 'Failed to generate random conformation',
1190 & ', itrial=',itrial
1191 write (*,*) 'Processor:',me,
1192 & ' Failed to generate random conformation',
1202 write (iout,'(a,i3,a)') 'Processor:',me,
1203 & ' error in generating random conformation.'
1204 write (*,'(a,i3,a)') 'Processor:',me,
1205 & ' error in generating random conformation.'
1208 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1213 & ' error in generating random conformation.'
1218 elseif (modecalc.eq.4) then
1219 read (inp,'(a)') intinname
1220 open (intin,file=intinname,status='old',err=333)
1221 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1222 & write (iout,'(a)') 'intinname',intinname
1223 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1225 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1227 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1229 stop 'Error opening angle file.'
1233 C Generate distance constraints, if the PDB structure is to be regularized.
1234 if (nthread.gt.0) then
1235 call read_threadbase
1238 if (me.eq.king .or. .not. out1file)
1240 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1241 write (iout,'(/a,i3,a)')
1242 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1243 write (iout,'(20i4)') (iss(i),i=1,ns)
1245 write(iout,*)"Running with dynamic disulfide-bond formation"
1247 write (iout,'(/a/)') 'Pre-formed links are:'
1253 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1254 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1260 if (ns.gt.0.and.dyn_ss) then
1264 forcon(i-nss)=forcon(i)
1271 dyn_ss_mask(iss(i))=.true.
1274 if (i2ndstr.gt.0) call secstrp2dihc
1275 c call geom_to_var(nvar,x)
1276 c call etotal(energia(0))
1277 c call enerprint(energia(0))
1278 c call briefout(0,etot)
1280 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1281 cd write (iout,'(a)') 'Variable list:'
1282 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1284 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1285 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1286 & 'Processor',myrank,': end reading molecular data.'
1290 c--------------------------------------------------------------------------
1291 logical function seq_comp(itypea,itypeb,length)
1293 integer length,itypea(length),itypeb(length)
1296 if (itypea(i).ne.itypeb(i)) then
1304 c-----------------------------------------------------------------------------
1305 subroutine read_bridge
1306 C Read information about disulfide bridges.
1307 implicit real*8 (a-h,o-z)
1308 include 'DIMENSIONS'
1312 include 'COMMON.IOUNITS'
1313 include 'COMMON.GEO'
1314 include 'COMMON.VAR'
1315 include 'COMMON.INTERACT'
1316 include 'COMMON.LOCAL'
1317 include 'COMMON.NAMES'
1318 include 'COMMON.CHAIN'
1319 include 'COMMON.FFIELD'
1320 include 'COMMON.SBRIDGE'
1321 include 'COMMON.HEADER'
1322 include 'COMMON.CONTROL'
1323 include 'COMMON.DBASE'
1324 include 'COMMON.THREAD'
1325 include 'COMMON.TIME1'
1326 include 'COMMON.SETUP'
1327 C Read bridging residues.
1328 read (inp,*) ns,(iss(i),i=1,ns)
1330 if(me.eq.king.or..not.out1file)
1331 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1332 C Check whether the specified bridging residues are cystines.
1334 if (itype(iss(i)).ne.1) then
1335 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1336 & 'Do you REALLY think that the residue ',
1337 & restyp(itype(iss(i))),i,
1338 & ' can form a disulfide bridge?!!!'
1339 write (*,'(2a,i3,a)')
1340 & 'Do you REALLY think that the residue ',
1341 & restyp(itype(iss(i))),i,
1342 & ' can form a disulfide bridge?!!!'
1344 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1349 C Read preformed bridges.
1351 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1353 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1356 C Check if the residues involved in bridges are in the specified list of
1357 C bridging residues.
1360 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1361 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1362 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1363 & ' contains residues present in other pairs.'
1364 write (*,'(a,i3,a)') 'Disulfide pair',i,
1365 & ' contains residues present in other pairs.'
1367 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1373 if (ihpb(i).eq.iss(j)) goto 10
1375 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1378 if (jhpb(i).eq.iss(j)) goto 20
1380 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1386 ihpb(i)=ihpb(i)+nres
1387 jhpb(i)=jhpb(i)+nres
1393 c----------------------------------------------------------------------------
1394 subroutine read_x(kanal,*)
1395 implicit real*8 (a-h,o-z)
1396 include 'DIMENSIONS'
1397 include 'COMMON.GEO'
1398 include 'COMMON.VAR'
1399 include 'COMMON.CHAIN'
1400 include 'COMMON.IOUNITS'
1401 include 'COMMON.CONTROL'
1402 include 'COMMON.LOCAL'
1403 include 'COMMON.INTERACT'
1404 c Read coordinates from input
1406 read(kanal,'(8f10.5)',end=10,err=10)
1407 & ((c(l,k),l=1,3),k=1,nres),
1408 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1411 c(j,2*nres)=c(j,nres)
1413 call int_from_cart1(.false.)
1416 dc(j,i)=c(j,i+1)-c(j,i)
1417 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1421 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1423 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1424 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1432 c----------------------------------------------------------------------------
1433 subroutine read_threadbase
1434 implicit real*8 (a-h,o-z)
1435 include 'DIMENSIONS'
1436 include 'COMMON.IOUNITS'
1437 include 'COMMON.GEO'
1438 include 'COMMON.VAR'
1439 include 'COMMON.INTERACT'
1440 include 'COMMON.LOCAL'
1441 include 'COMMON.NAMES'
1442 include 'COMMON.CHAIN'
1443 include 'COMMON.FFIELD'
1444 include 'COMMON.SBRIDGE'
1445 include 'COMMON.HEADER'
1446 include 'COMMON.CONTROL'
1447 include 'COMMON.DBASE'
1448 include 'COMMON.THREAD'
1449 include 'COMMON.TIME1'
1450 C Read pattern database for threading.
1451 read (icbase,*) nseq
1453 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1454 & nres_base(2,i),nres_base(3,i)
1455 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1457 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1458 c & nres_base(2,i),nres_base(3,i)
1459 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1463 if (weidis.eq.0.0D0) weidis=0.1D0
1472 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1473 write (iout,'(a,i5)') 'nexcl: ',nexcl
1474 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1477 c------------------------------------------------------------------------------
1478 subroutine setup_var
1479 implicit real*8 (a-h,o-z)
1480 include 'DIMENSIONS'
1481 include 'COMMON.IOUNITS'
1482 include 'COMMON.GEO'
1483 include 'COMMON.VAR'
1484 include 'COMMON.INTERACT'
1485 include 'COMMON.LOCAL'
1486 include 'COMMON.NAMES'
1487 include 'COMMON.CHAIN'
1488 include 'COMMON.FFIELD'
1489 include 'COMMON.SBRIDGE'
1490 include 'COMMON.HEADER'
1491 include 'COMMON.CONTROL'
1492 include 'COMMON.DBASE'
1493 include 'COMMON.THREAD'
1494 include 'COMMON.TIME1'
1495 C Set up variable list.
1501 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1503 ialph(i,1)=nvar+nside
1507 if (indphi.gt.0) then
1509 else if (indback.gt.0) then
1514 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1517 c----------------------------------------------------------------------------
1518 subroutine gen_dist_constr
1519 C Generate CA distance constraints.
1520 implicit real*8 (a-h,o-z)
1521 include 'DIMENSIONS'
1522 include 'COMMON.IOUNITS'
1523 include 'COMMON.GEO'
1524 include 'COMMON.VAR'
1525 include 'COMMON.INTERACT'
1526 include 'COMMON.LOCAL'
1527 include 'COMMON.NAMES'
1528 include 'COMMON.CHAIN'
1529 include 'COMMON.FFIELD'
1530 include 'COMMON.SBRIDGE'
1531 include 'COMMON.HEADER'
1532 include 'COMMON.CONTROL'
1533 include 'COMMON.DBASE'
1534 include 'COMMON.THREAD'
1535 include 'COMMON.TIME1'
1536 dimension itype_pdb(maxres)
1537 common /pizda/ itype_pdb
1539 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1540 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1541 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1543 do i=nstart_sup,nstart_sup+nsup-1
1544 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1545 cd & ' seq_pdb', restyp(itype_pdb(i))
1546 do j=i+2,nstart_sup+nsup-1
1548 ihpb(nhpb)=i+nstart_seq-nstart_sup
1549 jhpb(nhpb)=j+nstart_seq-nstart_sup
1551 dhpb(nhpb)=dist(i,j)
1554 cd write (iout,'(a)') 'Distance constraints:'
1559 cd if (ii.gt.nres) then
1564 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1565 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1566 cd & dhpb(i),forcon(i)
1570 c----------------------------------------------------------------------------
1572 implicit real*8 (a-h,o-z)
1573 include 'DIMENSIONS'
1574 include 'COMMON.MAP'
1575 include 'COMMON.IOUNITS'
1576 character*3 angid(4) /'THE','PHI','ALP','OME'/
1577 character*80 mapcard,ucase
1579 read (inp,'(a)') mapcard
1580 mapcard=ucase(mapcard)
1581 if (index(mapcard,'PHI').gt.0) then
1583 else if (index(mapcard,'THE').gt.0) then
1585 else if (index(mapcard,'ALP').gt.0) then
1587 else if (index(mapcard,'OME').gt.0) then
1590 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1591 stop 'Error - illegal variable spec in MAP card.'
1593 call readi (mapcard,'RES1',res1(imap),0)
1594 call readi (mapcard,'RES2',res2(imap),0)
1595 if (res1(imap).eq.0) then
1596 res1(imap)=res2(imap)
1597 else if (res2(imap).eq.0) then
1598 res2(imap)=res1(imap)
1600 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1602 & 'Error - illegal definition of variable group in MAP.'
1603 stop 'Error - illegal definition of variable group in MAP.'
1605 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1606 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1607 call readi(mapcard,'NSTEP',nstep(imap),0)
1608 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1610 & 'Illegal boundary and/or step size specification in MAP.'
1611 stop 'Illegal boundary and/or step size specification in MAP.'
1616 c----------------------------------------------------------------------------
1618 implicit real*8 (a-h,o-z)
1619 include 'DIMENSIONS'
1620 include 'COMMON.IOUNITS'
1621 include 'COMMON.GEO'
1622 include 'COMMON.CSA'
1623 include 'COMMON.BANK'
1624 include 'COMMON.CONTROL'
1626 character*620 mcmcard
1627 call card_concat(mcmcard)
1629 call readi(mcmcard,'NCONF',nconf,50)
1630 call readi(mcmcard,'NADD',nadd,0)
1631 call readi(mcmcard,'JSTART',jstart,1)
1632 call readi(mcmcard,'JEND',jend,1)
1633 call readi(mcmcard,'NSTMAX',nstmax,500000)
1634 call readi(mcmcard,'N0',n0,1)
1635 call readi(mcmcard,'N1',n1,6)
1636 call readi(mcmcard,'N2',n2,4)
1637 call readi(mcmcard,'N3',n3,0)
1638 call readi(mcmcard,'N4',n4,0)
1639 call readi(mcmcard,'N5',n5,0)
1640 call readi(mcmcard,'N6',n6,10)
1641 call readi(mcmcard,'N7',n7,0)
1642 call readi(mcmcard,'N8',n8,0)
1643 call readi(mcmcard,'N9',n9,0)
1644 call readi(mcmcard,'N14',n14,0)
1645 call readi(mcmcard,'N15',n15,0)
1646 call readi(mcmcard,'N16',n16,0)
1647 call readi(mcmcard,'N17',n17,0)
1648 call readi(mcmcard,'N18',n18,0)
1650 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1652 call readi(mcmcard,'NDIFF',ndiff,2)
1653 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1654 call readi(mcmcard,'IS1',is1,1)
1655 call readi(mcmcard,'IS2',is2,8)
1656 call readi(mcmcard,'NRAN0',nran0,4)
1657 call readi(mcmcard,'NRAN1',nran1,2)
1658 call readi(mcmcard,'IRR',irr,1)
1659 call readi(mcmcard,'NSEED',nseed,20)
1660 call readi(mcmcard,'NTOTAL',ntotal,10000)
1661 call reada(mcmcard,'CUT1',cut1,2.0d0)
1662 call reada(mcmcard,'CUT2',cut2,5.0d0)
1663 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1664 call readi(mcmcard,'ICMAX',icmax,3)
1665 call readi(mcmcard,'IRESTART',irestart,0)
1666 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1669 call reada(mcmcard,'DELE',dele,20.0d0)
1670 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1671 call readi(mcmcard,'IREF',iref,0)
1672 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1673 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1674 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1675 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1676 write (iout,*) "NCONF_IN",nconf_in
1679 c----------------------------------------------------------------------------
1680 cfmc subroutine mcmfread
1681 cfmc implicit real*8 (a-h,o-z)
1682 cfmc include 'DIMENSIONS'
1683 cfmc include 'COMMON.MCMF'
1684 cfmc include 'COMMON.IOUNITS'
1685 cfmc include 'COMMON.GEO'
1686 cfmc character*80 ucase
1687 cfmc character*620 mcmcard
1688 cfmc call card_concat(mcmcard)
1690 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1691 cfmc write(iout,*)'MAXRANT=',maxrant
1692 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1693 cfmc write(iout,*)'MAXFAM=',maxfam
1694 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1695 cfmc write(iout,*)'NNET1=',nnet1
1696 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1697 cfmc write(iout,*)'NNET2=',nnet2
1698 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1699 cfmc write(iout,*)'NNET3=',nnet3
1700 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1701 cfmc write(iout,*)'ILASTT=',ilastt
1702 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1703 cfmc write(iout,*)'MAXSTR=',maxstr
1704 cfmc maxstr_f=maxstr/maxfam
1705 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1706 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1707 cfmc write(iout,*)'NMCMF=',nmcmf
1708 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1709 cfmc write(iout,*)'IFOCUS=',ifocus
1710 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1711 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1712 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1713 cfmc write(iout,*)'INTPRT=',intprt
1714 cfmc call readi(mcmcard,'IPRT',iprt,100)
1715 cfmc write(iout,*)'IPRT=',iprt
1716 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1717 cfmc write(iout,*)'IMAXTR=',imaxtr
1718 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1719 cfmc write(iout,*)'MAXEVEN=',maxeven
1720 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1721 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1722 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1723 cfmc write(iout,*)'INIMIN=',inimin
1724 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1725 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1726 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1727 cfmc write(iout,*)'NTHREAD=',nthread
1728 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1729 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1730 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1731 cfmc write(iout,*)'MAXPERT=',maxpert
1732 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1733 cfmc write(iout,*)'IRMSD=',irmsd
1734 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1735 cfmc write(iout,*)'DENEMIN=',denemin
1736 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1737 cfmc write(iout,*)'RCUT1S=',rcut1s
1738 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1739 cfmc write(iout,*)'RCUT1E=',rcut1e
1740 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1741 cfmc write(iout,*)'RCUT2S=',rcut2s
1742 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1743 cfmc write(iout,*)'RCUT2E=',rcut2e
1744 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1745 cfmc write(iout,*)'DPERT1=',d_pert1
1746 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1747 cfmc write(iout,*)'DPERT1A=',d_pert1a
1748 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1749 cfmc write(iout,*)'DPERT2=',d_pert2
1750 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1751 cfmc write(iout,*)'DPERT2A=',d_pert2a
1752 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1753 cfmc write(iout,*)'DPERT2B=',d_pert2b
1754 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1755 cfmc write(iout,*)'DPERT2C=',d_pert2c
1756 cfmc d_pert1=deg2rad*d_pert1
1757 cfmc d_pert1a=deg2rad*d_pert1a
1758 cfmc d_pert2=deg2rad*d_pert2
1759 cfmc d_pert2a=deg2rad*d_pert2a
1760 cfmc d_pert2b=deg2rad*d_pert2b
1761 cfmc d_pert2c=deg2rad*d_pert2c
1762 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1763 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1764 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1765 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1766 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1767 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1768 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1769 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1770 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1771 cfmc write(iout,*)'RCUTINI=',rcutini
1772 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1773 cfmc write(iout,*)'GRAT=',grat
1774 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1775 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1779 c----------------------------------------------------------------------------
1781 implicit real*8 (a-h,o-z)
1782 include 'DIMENSIONS'
1783 include 'COMMON.MCM'
1784 include 'COMMON.MCE'
1785 include 'COMMON.IOUNITS'
1787 character*320 mcmcard
1788 call card_concat(mcmcard)
1789 call readi(mcmcard,'MAXACC',maxacc,100)
1790 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1791 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1792 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1793 call readi(mcmcard,'MAXREPM',maxrepm,200)
1794 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1795 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1796 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1797 call reada(mcmcard,'E_UP',e_up,5.0D0)
1798 call reada(mcmcard,'DELTE',delte,0.1D0)
1799 call readi(mcmcard,'NSWEEP',nsweep,5)
1800 call readi(mcmcard,'NSTEPH',nsteph,0)
1801 call readi(mcmcard,'NSTEPC',nstepc,0)
1802 call reada(mcmcard,'TMIN',tmin,298.0D0)
1803 call reada(mcmcard,'TMAX',tmax,298.0D0)
1804 call readi(mcmcard,'NWINDOW',nwindow,0)
1805 call readi(mcmcard,'PRINT_MC',print_mc,0)
1806 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1807 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1808 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1809 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1810 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1811 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1812 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1813 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1814 if (nwindow.gt.0) then
1815 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1817 winlen(i)=winend(i)-winstart(i)+1
1820 if (tmax.lt.tmin) tmax=tmin
1821 if (tmax.eq.tmin) then
1825 if (nstepc.gt.0 .and. nsteph.gt.0) then
1826 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1827 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1829 C Probabilities of different move types
1830 sumpro_type(0)=0.0D0
1831 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1832 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1833 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1834 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1835 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1836 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1837 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1839 print *,'i',i,' sumprotype',sumpro_type(i)
1840 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1841 print *,'i',i,' sumprotype',sumpro_type(i)
1845 c----------------------------------------------------------------------------
1846 subroutine read_minim
1847 implicit real*8 (a-h,o-z)
1848 include 'DIMENSIONS'
1849 include 'COMMON.MINIM'
1850 include 'COMMON.IOUNITS'
1851 include 'COMMON.CONTROL'
1852 include 'COMMON.SETUP'
1854 character*320 minimcard
1855 call card_concat(minimcard)
1856 call readi(minimcard,'MAXMIN',maxmin,2000)
1857 call readi(minimcard,'MAXFUN',maxfun,5000)
1858 call readi(minimcard,'MINMIN',minmin,maxmin)
1859 call readi(minimcard,'MINFUN',minfun,maxmin)
1860 call reada(minimcard,'TOLF',tolf,1.0D-2)
1861 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1862 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1863 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1864 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1866 if (.not. out1file .or. me.eq.king) then
1868 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1869 & 'Options in energy minimization:'
1870 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1871 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1872 & 'MinMin:',MinMin,' MinFun:',MinFun,
1873 & ' TolF:',TolF,' RTolF:',RTolF
1879 c----------------------------------------------------------------------------
1880 subroutine read_angles(kanal,*)
1881 implicit real*8 (a-h,o-z)
1882 include 'DIMENSIONS'
1883 include 'COMMON.GEO'
1884 include 'COMMON.VAR'
1885 include 'COMMON.CHAIN'
1886 include 'COMMON.IOUNITS'
1887 include 'COMMON.CONTROL'
1888 c Read angles from input
1890 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1891 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1892 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1893 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1896 c 9/7/01 avoid 180 deg valence angle
1897 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1899 theta(i)=deg2rad*theta(i)
1900 phi(i)=deg2rad*phi(i)
1901 alph(i)=deg2rad*alph(i)
1902 omeg(i)=deg2rad*omeg(i)
1907 c----------------------------------------------------------------------------
1908 subroutine reada(rekord,lancuch,wartosc,default)
1910 character*(*) rekord,lancuch
1911 double precision wartosc,default
1914 iread=index(rekord,lancuch)
1915 if (iread.eq.0) then
1919 iread=iread+ilen(lancuch)+1
1920 read (rekord(iread:),*,err=10,end=10) wartosc
1925 c----------------------------------------------------------------------------
1926 subroutine readi(rekord,lancuch,wartosc,default)
1928 character*(*) rekord,lancuch
1929 integer wartosc,default
1932 iread=index(rekord,lancuch)
1933 if (iread.eq.0) then
1937 iread=iread+ilen(lancuch)+1
1938 read (rekord(iread:),*,err=10,end=10) wartosc
1943 c----------------------------------------------------------------------------
1944 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1947 integer tablica(dim),default
1948 character*(*) rekord,lancuch
1955 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1956 if (iread.eq.0) return
1957 iread=iread+ilen(lancuch)+1
1958 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1961 c----------------------------------------------------------------------------
1962 subroutine multreada(rekord,lancuch,tablica,dim,default)
1965 double precision tablica(dim),default
1966 character*(*) rekord,lancuch
1973 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1974 if (iread.eq.0) return
1975 iread=iread+ilen(lancuch)+1
1976 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1979 c----------------------------------------------------------------------------
1980 subroutine openunits
1981 implicit real*8 (a-h,o-z)
1982 include 'DIMENSIONS'
1985 character*16 form,nodename
1988 include 'COMMON.SETUP'
1989 include 'COMMON.IOUNITS'
1991 include 'COMMON.CONTROL'
1992 integer lenpre,lenpot,ilen,lentmp
1994 character*3 out1file_text,ucase
1997 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1998 call getenv_loc("PREFIX",prefix)
2000 call getenv_loc("POT",pot)
2001 call getenv_loc("DIRTMP",tmpdir)
2002 call getenv_loc("CURDIR",curdir)
2003 call getenv_loc("OUT1FILE",out1file_text)
2004 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2005 out1file_text=ucase(out1file_text)
2006 if (out1file_text(1:1).eq."Y") then
2009 out1file=fg_rank.gt.0
2014 if (lentmp.gt.0) then
2015 write (*,'(80(1h!))')
2016 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2017 write (*,'(80(1h!))')
2018 write (*,*)"All output files will be on node /tmp directory."
2020 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2021 if (me.eq.king) then
2022 write (*,*) "The master node is ",nodename
2023 else if (fg_rank.eq.0) then
2024 write (*,*) "I am the CG slave node ",nodename
2026 write (*,*) "I am the FG slave node ",nodename
2029 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2030 lenpre = lentmp+lenpre+1
2032 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2033 C Get the names and open the input files
2034 #if defined(WINIFL) || defined(WINPGI)
2035 open(1,file=pref_orig(:ilen(pref_orig))//
2036 & '.inp',status='old',readonly,shared)
2037 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2038 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2039 C Get parameter filenames and open the parameter files.
2040 call getenv_loc('BONDPAR',bondname)
2041 open (ibond,file=bondname,status='old',readonly,shared)
2042 call getenv_loc('THETPAR',thetname)
2043 open (ithep,file=thetname,status='old',readonly,shared)
2045 call getenv_loc('THETPARPDB',thetname_pdb)
2046 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2048 call getenv_loc('ROTPAR',rotname)
2049 open (irotam,file=rotname,status='old',readonly,shared)
2051 call getenv_loc('ROTPARPDB',rotname_pdb)
2052 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2054 call getenv_loc('TORPAR',torname)
2055 open (itorp,file=torname,status='old',readonly,shared)
2056 call getenv_loc('TORDPAR',tordname)
2057 open (itordp,file=tordname,status='old',readonly,shared)
2058 call getenv_loc('FOURIER',fouriername)
2059 open (ifourier,file=fouriername,status='old',readonly,shared)
2060 call getenv_loc('ELEPAR',elename)
2061 open (ielep,file=elename,status='old',readonly,shared)
2062 call getenv_loc('SIDEPAR',sidename)
2063 open (isidep,file=sidename,status='old',readonly,shared)
2064 #elif (defined CRAY) || (defined AIX)
2065 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2067 c print *,"Processor",myrank," opened file 1"
2068 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2069 c print *,"Processor",myrank," opened file 9"
2070 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2071 C Get parameter filenames and open the parameter files.
2072 call getenv_loc('BONDPAR',bondname)
2073 open (ibond,file=bondname,status='old',action='read')
2074 c print *,"Processor",myrank," opened file IBOND"
2075 call getenv_loc('THETPAR',thetname)
2076 open (ithep,file=thetname,status='old',action='read')
2077 c print *,"Processor",myrank," opened file ITHEP"
2079 call getenv_loc('THETPARPDB',thetname_pdb)
2080 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2082 call getenv_loc('ROTPAR',rotname)
2083 open (irotam,file=rotname,status='old',action='read')
2084 c print *,"Processor",myrank," opened file IROTAM"
2086 call getenv_loc('ROTPARPDB',rotname_pdb)
2087 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2089 call getenv_loc('TORPAR',torname)
2090 open (itorp,file=torname,status='old',action='read')
2091 c print *,"Processor",myrank," opened file ITORP"
2092 call getenv_loc('TORDPAR',tordname)
2093 open (itordp,file=tordname,status='old',action='read')
2094 c print *,"Processor",myrank," opened file ITORDP"
2095 call getenv_loc('SCCORPAR',sccorname)
2096 open (isccor,file=sccorname,status='old',action='read')
2097 c print *,"Processor",myrank," opened file ISCCOR"
2098 call getenv_loc('FOURIER',fouriername)
2099 open (ifourier,file=fouriername,status='old',action='read')
2100 c print *,"Processor",myrank," opened file IFOURIER"
2101 call getenv_loc('ELEPAR',elename)
2102 open (ielep,file=elename,status='old',action='read')
2103 c print *,"Processor",myrank," opened file IELEP"
2104 call getenv_loc('SIDEPAR',sidename)
2105 open (isidep,file=sidename,status='old',action='read')
2106 c print *,"Processor",myrank," opened file ISIDEP"
2107 c print *,"Processor",myrank," opened parameter files"
2109 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2110 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2111 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2112 C Get parameter filenames and open the parameter files.
2113 call getenv_loc('BONDPAR',bondname)
2114 open (ibond,file=bondname,status='old')
2115 call getenv_loc('THETPAR',thetname)
2116 open (ithep,file=thetname,status='old')
2118 call getenv_loc('THETPARPDB',thetname_pdb)
2119 open (ithep_pdb,file=thetname_pdb,status='old')
2121 call getenv_loc('ROTPAR',rotname)
2122 open (irotam,file=rotname,status='old')
2124 call getenv_loc('ROTPARPDB',rotname_pdb)
2125 open (irotam_pdb,file=rotname_pdb,status='old')
2127 call getenv_loc('TORPAR',torname)
2128 open (itorp,file=torname,status='old')
2129 call getenv_loc('TORDPAR',tordname)
2130 open (itordp,file=tordname,status='old')
2131 call getenv_loc('SCCORPAR',sccorname)
2132 open (isccor,file=sccorname,status='old')
2133 call getenv_loc('FOURIER',fouriername)
2134 open (ifourier,file=fouriername,status='old')
2135 call getenv_loc('ELEPAR',elename)
2136 open (ielep,file=elename,status='old')
2137 call getenv_loc('SIDEPAR',sidename)
2138 open (isidep,file=sidename,status='old')
2139 call getenv_loc('LIPTRANPAR',liptranname)
2140 open (iliptranpar,file=liptranname,status='old')
2142 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2144 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2145 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2146 C Get parameter filenames and open the parameter files.
2147 call getenv_loc('BONDPAR',bondname)
2148 open (ibond,file=bondname,status='old',readonly)
2149 call getenv_loc('THETPAR',thetname)
2150 open (ithep,file=thetname,status='old',readonly)
2151 call getenv_loc('ROTPAR',rotname)
2152 open (irotam,file=rotname,status='old',readonly)
2153 call getenv_loc('TORPAR',torname)
2154 open (itorp,file=torname,status='old',readonly)
2155 call getenv_loc('TORDPAR',tordname)
2156 open (itordp,file=tordname,status='old',readonly)
2157 call getenv_loc('SCCORPAR',sccorname)
2158 open (isccor,file=sccorname,status='old',readonly)
2160 call getenv_loc('THETPARPDB',thetname_pdb)
2161 c print *,"thetname_pdb ",thetname_pdb
2162 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2163 c print *,ithep_pdb," opened"
2165 call getenv_loc('FOURIER',fouriername)
2166 open (ifourier,file=fouriername,status='old',readonly)
2167 call getenv_loc('ELEPAR',elename)
2168 open (ielep,file=elename,status='old',readonly)
2169 call getenv_loc('SIDEPAR',sidename)
2170 open (isidep,file=sidename,status='old',readonly)
2171 call getenv_loc('LIPTRANPAR',liptranname)
2172 open (iliptranpar,file=liptranname,status='old',readonly)
2174 call getenv_loc('ROTPARPDB',rotname_pdb)
2175 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2180 C 8/9/01 In the newest version SCp interaction constants are read from a file
2181 C Use -DOLDSCP to use hard-coded constants instead.
2183 call getenv_loc('SCPPAR',scpname)
2184 #if defined(WINIFL) || defined(WINPGI)
2185 open (iscpp,file=scpname,status='old',readonly,shared)
2186 #elif (defined CRAY) || (defined AIX)
2187 open (iscpp,file=scpname,status='old',action='read')
2189 open (iscpp,file=scpname,status='old')
2191 open (iscpp,file=scpname,status='old',readonly)
2194 call getenv_loc('PATTERN',patname)
2195 #if defined(WINIFL) || defined(WINPGI)
2196 open (icbase,file=patname,status='old',readonly,shared)
2197 #elif (defined CRAY) || (defined AIX)
2198 open (icbase,file=patname,status='old',action='read')
2200 open (icbase,file=patname,status='old')
2202 open (icbase,file=patname,status='old',readonly)
2205 C Open output file only for CG processes
2206 c print *,"Processor",myrank," fg_rank",fg_rank
2207 if (fg_rank.eq.0) then
2209 if (nodes.eq.1) then
2212 npos = dlog10(dfloat(nodes-1))+1
2214 if (npos.lt.3) npos=3
2215 write (liczba,'(i1)') npos
2216 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2218 write (liczba,form) me
2219 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2220 & liczba(:ilen(liczba))
2221 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2223 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2225 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2226 & liczba(:ilen(liczba))//'.mol2'
2227 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2228 & liczba(:ilen(liczba))//'.stat'
2230 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2231 & //liczba(:ilen(liczba))//'.stat')
2232 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2235 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2236 & liczba(:ilen(liczba))//'.const'
2241 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2242 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2243 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2244 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2245 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2247 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2249 rest2name=prefix(:ilen(prefix))//'.rst'
2251 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2254 #if defined(AIX) || defined(PGI)
2255 if (me.eq.king .or. .not. out1file)
2256 & open(iout,file=outname,status='unknown')
2258 if (fg_rank.gt.0) then
2259 write (liczba,'(i3.3)') myrank/nfgtasks
2260 write (ll,'(bz,i3.3)') fg_rank
2261 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2266 open(igeom,file=intname,status='unknown',position='append')
2267 open(ipdb,file=pdbname,status='unknown')
2268 open(imol2,file=mol2name,status='unknown')
2269 open(istat,file=statname,status='unknown',position='append')
2271 c1out open(iout,file=outname,status='unknown')
2274 if (me.eq.king .or. .not.out1file)
2275 & open(iout,file=outname,status='unknown')
2277 if (fg_rank.gt.0) then
2278 write (liczba,'(i3.3)') myrank/nfgtasks
2279 write (ll,'(bz,i3.3)') fg_rank
2280 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2285 open(igeom,file=intname,status='unknown',access='append')
2286 open(ipdb,file=pdbname,status='unknown')
2287 open(imol2,file=mol2name,status='unknown')
2288 open(istat,file=statname,status='unknown',access='append')
2290 c1out open(iout,file=outname,status='unknown')
2293 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2294 csa_seed=prefix(:lenpre)//'.CSA.seed'
2295 csa_history=prefix(:lenpre)//'.CSA.history'
2296 csa_bank=prefix(:lenpre)//'.CSA.bank'
2297 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2298 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2299 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2300 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2301 csa_int=prefix(:lenpre)//'.int'
2302 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2303 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2304 csa_in=prefix(:lenpre)//'.CSA.in'
2305 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2308 write (iout,'(80(1h-))')
2309 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2310 write (iout,'(80(1h-))')
2311 write (iout,*) "Input file : ",
2312 & pref_orig(:ilen(pref_orig))//'.inp'
2313 write (iout,*) "Output file : ",
2314 & outname(:ilen(outname))
2316 write (iout,*) "Sidechain potential file : ",
2317 & sidename(:ilen(sidename))
2319 write (iout,*) "SCp potential file : ",
2320 & scpname(:ilen(scpname))
2322 write (iout,*) "Electrostatic potential file : ",
2323 & elename(:ilen(elename))
2324 write (iout,*) "Cumulant coefficient file : ",
2325 & fouriername(:ilen(fouriername))
2326 write (iout,*) "Torsional parameter file : ",
2327 & torname(:ilen(torname))
2328 write (iout,*) "Double torsional parameter file : ",
2329 & tordname(:ilen(tordname))
2330 write (iout,*) "SCCOR parameter file : ",
2331 & sccorname(:ilen(sccorname))
2332 write (iout,*) "Bond & inertia constant file : ",
2333 & bondname(:ilen(bondname))
2334 write (iout,*) "Bending parameter file : ",
2335 & thetname(:ilen(thetname))
2336 write (iout,*) "Rotamer parameter file : ",
2337 & rotname(:ilen(rotname))
2338 write (iout,*) "Thetpdb parameter file : ",
2339 & thetname_pdb(:ilen(thetname_pdb))
2340 write (iout,*) "Threading database : ",
2341 & patname(:ilen(patname))
2343 &write (iout,*)" DIRTMP : ",
2345 write (iout,'(80(1h-))')
2349 c----------------------------------------------------------------------------
2350 subroutine card_concat(card)
2351 implicit real*8 (a-h,o-z)
2352 include 'DIMENSIONS'
2353 include 'COMMON.IOUNITS'
2355 character*80 karta,ucase
2357 read (inp,'(a)') karta
2360 do while (karta(80:80).eq.'&')
2361 card=card(:ilen(card)+1)//karta(:79)
2362 read (inp,'(a)') karta
2365 card=card(:ilen(card)+1)//karta
2368 c----------------------------------------------------------------------------------
2370 implicit real*8 (a-h,o-z)
2371 include 'DIMENSIONS'
2372 include 'COMMON.CHAIN'
2373 include 'COMMON.IOUNITS'
2375 include 'COMMON.CONTROL'
2376 open(irest2,file=rest2name,status='unknown')
2377 read(irest2,*) totT,EK,potE,totE,t_bath
2380 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2383 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2385 if(usampl.or.homol_nset.gt.1) then
2386 read (irest2,*) iset
2391 c---------------------------------------------------------------------------------
2392 subroutine read_fragments
2393 implicit real*8 (a-h,o-z)
2394 include 'DIMENSIONS'
2398 include 'COMMON.SETUP'
2399 include 'COMMON.CHAIN'
2400 include 'COMMON.IOUNITS'
2402 include 'COMMON.CONTROL'
2404 read(inp,*) nset,nfrag,npair,nfrag_back
2405 if(me.eq.king.or..not.out1file)
2406 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2407 & " nfrag_back",nfrag_back
2409 read(inp,*) mset(iset1)
2411 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2413 if(me.eq.king.or..not.out1file)
2414 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2415 & ifrag(2,i,iset1), qinfrag(i,iset1)
2418 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2420 if(me.eq.king.or..not.out1file)
2421 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2422 & ipair(2,i,iset1), qinpair(i,iset1)
2425 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2426 & wfrag_back(3,i,iset1),
2427 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2428 if(me.eq.king.or..not.out1file)
2429 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2430 & wfrag_back(2,i,iset1),
2431 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2432 & ifrag_back(2,i,iset1)
2437 C---------------------------------------------------------------------------
2438 subroutine read_afminp
2439 implicit real*8 (a-h,o-z)
2440 include 'DIMENSIONS'
2444 include 'COMMON.SETUP'
2445 include 'COMMON.CONTROL'
2446 include 'COMMON.CHAIN'
2447 include 'COMMON.IOUNITS'
2448 include 'COMMON.SBRIDGE'
2449 character*320 afmcard
2450 c print *, "wchodze"
2451 call card_concat(afmcard)
2452 call readi(afmcard,"BEG",afmbeg,0)
2453 call readi(afmcard,"END",afmend,0)
2454 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2455 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2456 print *,'FORCE=' ,forceAFMconst
2457 CCCC NOW PROPERTIES FOR AFM
2460 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2462 distafminit=dsqrt(distafminit)
2463 print *,'initdist',distafminit
2467 c-------------------------------------------------------------------------------
2468 subroutine read_saxs_constr
2469 implicit real*8 (a-h,o-z)
2470 include 'DIMENSIONS'
2474 include 'COMMON.SETUP'
2475 include 'COMMON.CONTROL'
2476 include 'COMMON.CHAIN'
2477 include 'COMMON.IOUNITS'
2478 include 'COMMON.SBRIDGE'
2479 double precision cm(3)
2481 write (iout,*) "Calling read_saxs nsaxs",nsaxs
2483 if (saxs_mode.eq.0) then
2484 c SAXS distance distribution
2486 read(inp,*) distsaxs(i),Psaxs(i)
2490 Cnorm = Cnorm + Psaxs(i)
2492 write (iout,*) "Cnorm",Cnorm
2494 Psaxs(i)=Psaxs(i)/Cnorm
2496 write (iout,*) "Normalized distance distribution from SAXS"
2498 write (iout,'(f8.2,e15.5)') distsaxs(i),Psaxs(i)
2503 read (inp,'(30x,3f8.3)') (Csaxs(j,i),j=1,3)
2510 cm(j)=cm(j)+Csaxs(j,i)
2518 Csaxs(j,i)=Csaxs(j,i)-cm(j)
2521 write (iout,*) "SAXS sphere coordinates"
2523 write (iout,'(i5,3f10.5)') i,(Csaxs(j,i),j=1,3)
2528 c-------------------------------------------------------------------------------
2529 subroutine read_dist_constr
2530 implicit real*8 (a-h,o-z)
2531 include 'DIMENSIONS'
2535 include 'COMMON.SETUP'
2536 include 'COMMON.CONTROL'
2537 include 'COMMON.CHAIN'
2538 include 'COMMON.IOUNITS'
2539 include 'COMMON.SBRIDGE'
2540 integer ifrag_(2,100),ipair_(2,100)
2541 double precision wfrag_(100),wpair_(100)
2542 character*500 controlcard
2543 c write (iout,*) "Calling read_dist_constr"
2544 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2546 call card_concat(controlcard)
2547 call readi(controlcard,"NFRAG",nfrag_,0)
2548 call readi(controlcard,"NPAIR",npair_,0)
2549 call readi(controlcard,"NDIST",ndist_,0)
2550 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2551 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2552 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2553 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2554 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2555 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2556 c write (iout,*) "IFRAG"
2558 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2560 c write (iout,*) "IPAIR"
2562 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2564 if (.not.refstr .and. nfrag.gt.0) then
2566 & "ERROR: no reference structure to compute distance restraints"
2568 & "Restraints must be specified explicitly (NDIST=number)"
2571 if (nfrag.lt.2 .and. npair.gt.0) then
2572 write (iout,*) "ERROR: Less than 2 fragments specified",
2573 & " but distance restraints between pairs requested"
2578 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2579 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2580 & ifrag_(2,i)=nstart_sup+nsup-1
2581 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2583 if (wfrag_(i).gt.0.0d0) then
2584 do j=ifrag_(1,i),ifrag_(2,i)-1
2585 do k=j+1,ifrag_(2,i)
2586 c write (iout,*) "j",j," k",k
2588 if (constr_dist.eq.1) then
2593 forcon(nhpb)=wfrag_(i)
2594 else if (constr_dist.eq.2) then
2595 if (ddjk.le.dist_cut) then
2600 forcon(nhpb)=wfrag_(i)
2607 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2610 if (.not.out1file .or. me.eq.king)
2611 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2612 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2614 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2615 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2622 if (wpair_(i).gt.0.0d0) then
2630 do j=ifrag_(1,ii),ifrag_(2,ii)
2631 do k=ifrag_(1,jj),ifrag_(2,jj)
2635 forcon(nhpb)=wpair_(i)
2636 dhpb(nhpb)=dist(j,k)
2638 if (.not.out1file .or. me.eq.king)
2639 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2640 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2642 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2643 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2650 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2651 if (forcon(nhpb+1).gt.0.0d0) then
2653 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2655 if (.not.out1file .or. me.eq.king)
2656 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2657 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2659 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2660 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2667 c-------------------------------------------------------------------------------
2669 subroutine read_constr_homology
2671 include 'DIMENSIONS'
2675 include 'COMMON.SETUP'
2676 include 'COMMON.CONTROL'
2677 include 'COMMON.CHAIN'
2678 include 'COMMON.IOUNITS'
2680 include 'COMMON.GEO'
2681 include 'COMMON.INTERACT'
2682 include 'COMMON.NAMES'
2684 c For new homol impl
2686 include 'COMMON.VAR'
2689 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2691 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2692 c & sigma_odl_temp(maxres,maxres,max_template)
2694 character*24 model_ki_dist, model_ki_angle
2695 character*500 controlcard
2696 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2697 logical lprn /.true./
2702 c FP - Nov. 2014 Temporary specifications for new vars
2704 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2706 double precision, dimension (max_template,maxres) :: rescore
2707 double precision, dimension (max_template,maxres) :: rescore2
2708 double precision, dimension (max_template,maxres) :: rescore3
2709 character*24 pdbfile,tpl_k_rescore
2710 c -----------------------------------------------------------------
2711 c Reading multiple PDB ref structures and calculation of retraints
2712 c not using pre-computed ones stored in files model_ki_{dist,angle}
2714 c -----------------------------------------------------------------
2717 c Alternative: reading from input
2718 call card_concat(controlcard)
2719 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2720 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2721 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2722 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2723 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2724 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2725 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2726 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2727 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2728 if(.not.read2sigma.and.start_from_model) then
2729 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2730 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2731 start_from_model=.false.
2733 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2734 & write(iout,*) 'START_FROM_MODELS is ON'
2735 if(start_from_model .and. rest) then
2736 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2737 write(iout,*) 'START_FROM_MODELS is OFF'
2738 write(iout,*) 'remove restart keyword from input'
2741 if (homol_nset.gt.1)then
2742 call card_concat(controlcard)
2743 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2744 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2745 write(iout,*) "iset homology_weight "
2747 write(iout,*) i,waga_homology(i)
2750 iset=mod(kolor,homol_nset)+1
2753 waga_homology(1)=1.0
2756 cd write (iout,*) "nnt",nnt," nct",nct
2763 c write(iout,*) 'nnt=',nnt,'nct=',nct
2766 do k=1,constr_homology
2779 do k=1,constr_homology
2781 read(inp,'(a)') pdbfile
2782 if(me.eq.king .or. .not. out1file)
2783 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2784 & pdbfile(:ilen(pdbfile))
2785 open(ipdbin,file=pdbfile,status='old',err=33)
2787 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2788 & pdbfile(:ilen(pdbfile))
2791 c print *,'Begin reading pdb data'
2793 c Files containing res sim or local scores (former containing sigmas)
2796 write(kic2,'(bz,i2.2)') k
2798 tpl_k_rescore="template"//kic2//".sco"
2801 if (read2sigma) then
2802 call readpdb_template(k)
2807 c Distance restraints
2810 C Copy the coordinates from reference coordinates (?)
2814 c write (iout,*) "c(",j,i,") =",c(j,i)
2818 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2820 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2821 open (ientin,file=tpl_k_rescore,status='old')
2822 if (nnt.gt.1) rescore(k,1)=0.0d0
2823 do irec=nnt,nct ! loop for reading res sim
2824 if (read2sigma) then
2825 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2826 & rescore3_tmp,idomain_tmp
2828 idomain(k,i_tmp)=idomain_tmp
2829 rescore(k,i_tmp)=rescore_tmp
2830 rescore2(k,i_tmp)=rescore2_tmp
2831 rescore3(k,i_tmp)=rescore3_tmp
2832 if (.not. out1file .or. me.eq.king)
2833 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
2834 & i_tmp,rescore2_tmp,rescore_tmp,
2835 & rescore3_tmp,idomain_tmp
2838 read (ientin,*,end=1401) rescore_tmp
2840 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2841 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2842 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2847 if (waga_dist.ne.0.0d0) then
2855 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2856 c write (iout,*) k,i,j,distal,dist2_cut
2858 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2859 & .and. distal.le.dist2_cut ) then
2865 c write (iout,*) "k",k
2866 c write (iout,*) "i",i," j",j," constr_homology",
2871 if (read2sigma) then
2874 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2876 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2877 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2878 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2880 if (odl(k,ii).le.dist_cut) then
2881 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2884 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2885 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2887 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2888 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2892 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2895 l_homo(k,ii)=.false.
2902 c Theta, dihedral and SC retraints
2904 if (waga_angle.gt.0.0d0) then
2905 c open (ientin,file=tpl_k_sigma_dih,status='old')
2906 c do irec=1,maxres-3 ! loop for reading sigma_dih
2907 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2908 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2909 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2910 c & sigma_dih(k,i+nnt-1)
2915 if (idomain(k,i).eq.0) then
2919 dih(k,i)=phiref(i) ! right?
2920 c read (ientin,*) sigma_dih(k,i) ! original variant
2921 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2922 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2923 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2924 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2925 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2927 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2928 & rescore(k,i-2)+rescore(k,i-3))/4.0
2929 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2930 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2931 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2932 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2933 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2934 c Instead of res sim other local measure of b/b str reliability possible
2935 if (sigma_dih(k,i).ne.0)
2936 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2937 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2942 if (waga_theta.gt.0.0d0) then
2943 c open (ientin,file=tpl_k_sigma_theta,status='old')
2944 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2945 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2946 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2947 c & sigma_theta(k,i+nnt-1)
2952 do i = nnt+2,nct ! right? without parallel.
2953 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2954 c do i=ithet_start,ithet_end ! with FG parallel.
2955 if (idomain(k,i).eq.0) then
2956 sigma_theta(k,i)=0.0
2959 thetatpl(k,i)=thetaref(i)
2960 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2961 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2962 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2963 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2964 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2965 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2966 & rescore(k,i-2))/3.0
2967 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2968 if (sigma_theta(k,i).ne.0)
2969 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2971 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2972 c rescore(k,i-2) ! right expression ?
2973 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2977 if (waga_d.gt.0.0d0) then
2978 c open (ientin,file=tpl_k_sigma_d,status='old')
2979 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2980 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2981 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2982 c & sigma_d(k,i+nnt-1)
2986 do i = nnt,nct ! right? without parallel.
2987 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2988 c do i=loc_start,loc_end ! with FG parallel.
2989 if (itype(i).eq.10) cycle
2990 if (idomain(k,i).eq.0 ) then
2997 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2998 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2999 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3000 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
3001 sigma_d(k,i)=rescore3(k,i) ! right expression ?
3002 if (sigma_d(k,i).ne.0)
3003 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3005 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
3006 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3007 c read (ientin,*) sigma_d(k,i) ! 1st variant
3012 c remove distance restraints not used in any model from the list
3013 c shift data in all arrays
3015 if (waga_dist.ne.0.0d0) then
3021 if (ii_in_use(ii).eq.0.and.liiflag) then
3025 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
3026 & .not.liiflag.and.ii.eq.lim_odl) then
3027 if (ii.eq.lim_odl) then
3028 iishift=ii-iistart+1
3033 do ki=iistart,lim_odl-iishift
3034 ires_homo(ki)=ires_homo(ki+iishift)
3035 jres_homo(ki)=jres_homo(ki+iishift)
3036 ii_in_use(ki)=ii_in_use(ki+iishift)
3037 do k=1,constr_homology
3038 odl(k,ki)=odl(k,ki+iishift)
3039 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
3040 l_homo(k,ki)=l_homo(k,ki+iishift)
3044 lim_odl=lim_odl-iishift
3049 if (constr_homology.gt.0) call homology_partition
3050 if (constr_homology.gt.0) call init_int_table
3051 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3052 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3056 if (.not.lprn) return
3057 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3058 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3059 write (iout,*) "Distance restraints from templates"
3061 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
3062 & ii,ires_homo(ii),jres_homo(ii),
3063 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3064 & ki=1,constr_homology)
3066 write (iout,*) "Dihedral angle restraints from templates"
3068 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3069 & (rad2deg*dih(ki,i),
3070 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3072 write (iout,*) "Virtual-bond angle restraints from templates"
3074 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
3075 & (rad2deg*thetatpl(ki,i),
3076 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3078 write (iout,*) "SC restraints from templates"
3080 write(iout,'(i5,100(4f8.2,4x))') i,
3081 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3082 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3085 c -----------------------------------------------------------------
3088 c----------------------------------------------------------------------
3091 subroutine flush(iu)
3096 subroutine flush(iu)
3101 c------------------------------------------------------------------------------
3102 subroutine copy_to_tmp(source)
3103 include "DIMENSIONS"
3104 include "COMMON.IOUNITS"
3105 character*(*) source
3106 character* 256 tmpfile
3110 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3111 inquire(file=tmpfile,exist=ex)
3113 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3114 & " to temporary directory..."
3115 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3116 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3120 c------------------------------------------------------------------------------
3121 subroutine move_from_tmp(source)
3122 include "DIMENSIONS"
3123 include "COMMON.IOUNITS"
3124 character*(*) source
3127 write (*,*) "Moving ",source(:ilen(source)),
3128 & " from temporary directory to working directory"
3129 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3130 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3133 c------------------------------------------------------------------------------
3134 subroutine random_init(seed)
3136 C Initialize random number generator
3138 implicit real*8 (a-h,o-z)
3139 include 'DIMENSIONS'
3145 logical OKRandom, prng_restart
3147 integer iseed_array(4)
3149 include 'COMMON.IOUNITS'
3150 include 'COMMON.TIME1'
3151 include 'COMMON.THREAD'
3152 include 'COMMON.SBRIDGE'
3153 include 'COMMON.CONTROL'
3154 include 'COMMON.MCM'
3155 include 'COMMON.MAP'
3156 include 'COMMON.HEADER'
3157 include 'COMMON.CSA'
3158 include 'COMMON.CHAIN'
3159 include 'COMMON.MUCA'
3161 include 'COMMON.FFIELD'
3162 include 'COMMON.SETUP'
3163 iseed=-dint(dabs(seed))
3164 if (iseed.eq.0) then
3165 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3166 & 'Random seed undefined. The program will stop.'
3167 write (*,'(/80(1h*)/20x,a/80(1h*))')
3168 & 'Random seed undefined. The program will stop.'
3170 call mpi_finalize(mpi_comm_world,ierr)
3172 stop 'Bad random seed.'
3175 if (fg_rank.eq.0) then
3179 if(me.eq.king .or. .not. out1file)
3180 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3181 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3182 OKRandom = prng_restart(me,iseedi8)
3185 tmp=65536.0d0**(4-i)
3186 iseed_array(i) = dint(seed/tmp)
3187 seed=seed-iseed_array(i)*tmp
3190 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3191 & (iseed_array(i),i=1,4)
3192 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3193 & (iseed_array(i),i=1,4)
3194 OKRandom = prng_restart(me,iseed_array)
3197 c r1 = prng_next(me)
3198 r1=ran_number(0.0D0,1.0D0)
3199 if(me.eq.king .or. .not. out1file)
3200 & write (iout,*) 'ran_num',r1
3201 if (r1.lt.0.0d0) OKRandom=.false.
3203 if (.not.OKRandom) then
3204 write (iout,*) 'PRNG IS NOT WORKING!!!'
3205 print *,'PRNG IS NOT WORKING!!!'
3208 call mpi_abort(mpi_comm_world,error_msg,ierr)
3211 write (iout,*) 'too many processors for parallel prng'
3212 write (*,*) 'too many processors for parallel prng'
3220 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)