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,'CONSTR_HOMOL',constr_homology,0)
102 call readi(controlcard,'SYM',symetr,1)
103 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
104 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
105 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
106 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
107 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
108 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
109 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
110 call reada(controlcard,'DRMS',drms,0.1D0)
111 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
112 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
113 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
114 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
115 write (iout,'(a,f10.1)')'DRMS = ',drms
116 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
117 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
119 call readi(controlcard,'NZ_START',nz_start,0)
120 call readi(controlcard,'NZ_END',nz_end,0)
121 call readi(controlcard,'IZ_SC',iz_sc,0)
123 safety = 60.0d0*safety
126 call reada(controlcard,"T_BATH",t_bath,300.0d0)
127 minim=(index(controlcard,'MINIMIZE').gt.0)
128 dccart=(index(controlcard,'CART').gt.0)
129 overlapsc=(index(controlcard,'OVERLAP').gt.0)
130 overlapsc=.not.overlapsc
131 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
132 searchsc=.not.searchsc
133 sideadd=(index(controlcard,'SIDEADD').gt.0)
134 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
135 outpdb=(index(controlcard,'PDBOUT').gt.0)
136 outmol2=(index(controlcard,'MOL2OUT').gt.0)
137 pdbref=(index(controlcard,'PDBREF').gt.0)
138 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
139 indpdb=index(controlcard,'PDBSTART')
140 extconf=(index(controlcard,'EXTCONF').gt.0)
141 AFMlog=(index(controlcard,'AFM'))
142 selfguide=(index(controlcard,'SELFGUIDE'))
143 c print *,'AFMlog',AFMlog,selfguide,"KUPA"
144 call readi(controlcard,'IPRINT',iprint,0)
145 call readi(controlcard,'MAXGEN',maxgen,10000)
146 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
147 call readi(controlcard,"KDIAG",kdiag,0)
148 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
149 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
150 & write (iout,*) "RESCALE_MODE",rescale_mode
151 split_ene=index(controlcard,'SPLIT_ENE').gt.0
152 if (index(controlcard,'REGULAR').gt.0.0D0) then
153 call reada(controlcard,'WEIDIS',weidis,0.1D0)
157 if (index(controlcard,'CHECKGRAD').gt.0) then
159 if (index(controlcard,'CART').gt.0) then
161 elseif (index(controlcard,'CARINT').gt.0) then
166 elseif (index(controlcard,'THREAD').gt.0) then
168 call readi(controlcard,'THREAD',nthread,0)
169 if (nthread.gt.0) then
170 call reada(controlcard,'WEIDIS',weidis,0.1D0)
173 & write (iout,'(a)')'A number has to follow the THREAD keyword.'
174 stop 'Error termination in Read_Control.'
176 else if (index(controlcard,'MCMA').gt.0) then
178 else if (index(controlcard,'MCEE').gt.0) then
180 else if (index(controlcard,'MULTCONF').gt.0) then
182 else if (index(controlcard,'MAP').gt.0) then
184 call readi(controlcard,'MAP',nmap,0)
185 else if (index(controlcard,'CSA').gt.0) then
187 crc else if (index(controlcard,'ZSCORE').gt.0) then
189 crc ZSCORE is rm from UNRES, modecalc=9 is available
192 cfcm else if (index(controlcard,'MCMF').gt.0) then
194 else if (index(controlcard,'SOFTREG').gt.0) then
196 else if (index(controlcard,'CHECK_BOND').gt.0) then
198 else if (index(controlcard,'TEST').gt.0) then
200 else if (index(controlcard,'MD').gt.0) then
202 else if (index(controlcard,'RE ').gt.0) then
206 lmuca=index(controlcard,'MUCA').gt.0
207 call readi(controlcard,'MUCADYN',mucadyn,0)
208 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
209 if (lmuca .and. (me.eq.king .or. .not.out1file ))
211 write (iout,*) 'MUCADYN=',mucadyn
212 write (iout,*) 'MUCASMOOTH=',muca_smooth
215 iscode=index(controlcard,'ONE_LETTER')
216 indphi=index(controlcard,'PHI')
217 indback=index(controlcard,'BACK')
218 iranconf=index(controlcard,'RAND_CONF')
219 i2ndstr=index(controlcard,'USE_SEC_PRED')
220 gradout=index(controlcard,'GRADOUT').gt.0
221 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
222 C DISTCHAINMAX become obsolete for periodic boundry condition
223 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
224 C Reading the dimensions of box in x,y,z coordinates
225 call reada(controlcard,'BOXX',boxxsize,100.0d0)
226 call reada(controlcard,'BOXY',boxysize,100.0d0)
227 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
228 c Cutoff range for interactions
229 call reada(controlcard,"R_CUT",r_cut,15.0d0)
230 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
231 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
232 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
233 if (lipthick.gt.0.0d0) then
234 bordliptop=(boxzsize+lipthick)/2.0
235 bordlipbot=bordliptop-lipthick
237 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0))
238 & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
239 buflipbot=bordlipbot+lipbufthick
240 bufliptop=bordliptop-lipbufthick
241 if ((lipbufthick*2.0d0).gt.lipthick)
242 &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
244 c write(iout,*) "bordliptop=",bordliptop
245 c write(iout,*) "bordlipbot=",bordlipbot
246 c write(iout,*) "bufliptop=",bufliptop
247 c write(iout,*) "buflipbot=",buflipbot
250 if (me.eq.king .or. .not.out1file )
251 & write (iout,*) "DISTCHAINMAX",distchainmax
253 if(me.eq.king.or..not.out1file)
254 & write (iout,'(2a)') diagmeth(kdiag),
255 & ' routine used to diagonalize matrices.'
258 c--------------------------------------------------------------------------
259 subroutine read_REMDpar
263 implicit real*8 (a-h,o-z)
265 include 'COMMON.IOUNITS'
266 include 'COMMON.TIME1'
269 include 'COMMON.LANGEVIN'
271 include 'COMMON.LANGEVIN.lang0'
273 include 'COMMON.INTERACT'
274 include 'COMMON.NAMES'
276 include 'COMMON.REMD'
277 include 'COMMON.CONTROL'
278 include 'COMMON.SETUP'
280 character*320 controlcard
281 character*3200 controlcard1
282 integer iremd_m_total
284 if(me.eq.king.or..not.out1file)
285 & write (iout,*) "REMD setup"
287 call card_concat(controlcard)
288 call readi(controlcard,"NREP",nrep,3)
289 call readi(controlcard,"NSTEX",nstex,1000)
290 call reada(controlcard,"RETMIN",retmin,10.0d0)
291 call reada(controlcard,"RETMAX",retmax,1000.0d0)
292 mremdsync=(index(controlcard,'SYNC').gt.0)
293 call readi(controlcard,"NSYN",i_sync_step,100)
294 restart1file=(index(controlcard,'REST1FILE').gt.0)
295 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
296 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
297 if(max_cache_traj_use.gt.max_cache_traj)
298 & max_cache_traj_use=max_cache_traj
299 if(me.eq.king.or..not.out1file) then
300 cd if (traj1file) then
301 crc caching is in testing - NTWX is not ignored
302 cd write (iout,*) "NTWX value is ignored"
303 cd write (iout,*) " trajectory is stored to one file by master"
304 cd write (iout,*) " before exchange at NSTEX intervals"
306 write (iout,*) "NREP= ",nrep
307 write (iout,*) "NSTEX= ",nstex
308 write (iout,*) "SYNC= ",mremdsync
309 write (iout,*) "NSYN= ",i_sync_step
310 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
313 if (index(controlcard,'TLIST').gt.0) then
315 call card_concat(controlcard1)
316 read(controlcard1,*) (remd_t(i),i=1,nrep)
317 if(me.eq.king.or..not.out1file)
318 & write (iout,*)'tlist',(remd_t(i),i=1,nrep)
321 if (index(controlcard,'MLIST').gt.0) then
323 call card_concat(controlcard1)
324 read(controlcard1,*) (remd_m(i),i=1,nrep)
325 if(me.eq.king.or..not.out1file) then
326 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
329 iremd_m_total=iremd_m_total+remd_m(i)
331 write (iout,*) 'Total number of replicas ',iremd_m_total
334 if(me.eq.king.or..not.out1file)
335 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
338 c--------------------------------------------------------------------------
339 subroutine read_MDpar
343 implicit real*8 (a-h,o-z)
345 include 'COMMON.IOUNITS'
346 include 'COMMON.TIME1'
349 include 'COMMON.LANGEVIN'
351 include 'COMMON.LANGEVIN.lang0'
353 include 'COMMON.INTERACT'
354 include 'COMMON.NAMES'
356 include 'COMMON.SETUP'
357 include 'COMMON.CONTROL'
358 include 'COMMON.SPLITELE'
360 character*320 controlcard
362 call card_concat(controlcard)
363 call readi(controlcard,"NSTEP",n_timestep,1000000)
364 call readi(controlcard,"NTWE",ntwe,100)
365 call readi(controlcard,"NTWX",ntwx,1000)
366 call reada(controlcard,"DT",d_time,1.0d-1)
367 call reada(controlcard,"DVMAX",dvmax,2.0d1)
368 call reada(controlcard,"DAMAX",damax,1.0d1)
369 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
370 call readi(controlcard,"LANG",lang,0)
371 RESPA = index(controlcard,"RESPA") .gt. 0
372 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
373 ntime_split0=ntime_split
374 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
375 ntime_split0=ntime_split
376 c call reada(controlcard,"R_CUT",r_cut,2.0d0)
377 c call reada(controlcard,"LAMBDA",rlamb,0.3d0)
378 rest = index(controlcard,"REST").gt.0
379 tbf = index(controlcard,"TBF").gt.0
380 usampl = index(controlcard,"USAMPL").gt.0
382 mdpdb = index(controlcard,"MDPDB").gt.0
383 call reada(controlcard,"T_BATH",t_bath,300.0d0)
384 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
385 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
386 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
387 if (count_reset_moment.eq.0) count_reset_moment=1000000000
388 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
389 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
390 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
391 if (count_reset_vel.eq.0) count_reset_vel=1000000000
392 large = index(controlcard,"LARGE").gt.0
393 print_compon = index(controlcard,"PRINT_COMPON").gt.0
394 rattle = index(controlcard,"RATTLE").gt.0
395 preminim = index(controlcard,"PREMINIM").gt.0
397 dccart=(index(controlcard,'CART').gt.0)
400 c if performing umbrella sampling, fragments constrained are read from the fragment file
406 if(me.eq.king.or..not.out1file) then
408 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
410 write (iout,'(a)') "The units are:"
411 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
412 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
413 & " acceleration: angstrom/(48.9 fs)**2"
414 write (iout,'(a)') "energy: kcal/mol, temperature: K"
416 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
417 write (iout,'(a60,f10.5,a)')
418 & "Initial time step of numerical integration:",d_time,
420 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
422 write (iout,'(2a,i4,a)')
423 & "A-MTS algorithm used; initial time step for fast-varying",
424 & " short-range forces split into",ntime_split," steps."
425 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
426 & r_cut," lambda",rlamb
428 write (iout,'(2a,f10.5)')
429 & "Maximum acceleration threshold to reduce the time step",
430 & "/increase split number:",damax
431 write (iout,'(2a,f10.5)')
432 & "Maximum predicted energy drift to reduce the timestep",
433 & "/increase split number:",edriftmax
434 write (iout,'(a60,f10.5)')
435 & "Maximum velocity threshold to reduce velocities:",dvmax
436 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
437 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
438 if (rattle) write (iout,'(a60)')
439 & "Rattle algorithm used to constrain the virtual bonds"
440 if (preminim .or. iranconf.gt.0) then
442 & "Initial structure will be energy-minimized"
447 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
448 call reada(controlcard,"RWAT",rwat,1.4d0)
449 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
450 surfarea=index(controlcard,"SURFAREA").gt.0
451 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
452 if(me.eq.king.or..not.out1file)then
453 write (iout,'(/a,$)') "Langevin dynamics calculation"
456 & " with direct integration of Langevin equations"
457 else if (lang.eq.2) then
458 write (iout,'(a/)') " with TINKER stochasic MD integrator"
459 else if (lang.eq.3) then
460 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
461 else if (lang.eq.4) then
462 write (iout,'(a/)') " in overdamped mode"
464 write (iout,'(//a,i5)')
465 & "=========== ERROR: Unknown Langevin dynamics mode:",lang
468 write (iout,'(a60,f10.5)') "Temperature:",t_bath
469 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
470 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
471 write (iout,'(a60,f10.5)')
472 & "Scaling factor of the friction forces:",scal_fric
473 if (surfarea) write (iout,'(2a,i10,a)')
474 & "Friction coefficients will be scaled by solvent-accessible",
475 & " surface area every",reset_fricmat," steps."
477 c Calculate friction coefficients and bounds of stochastic forces
478 eta=6*pi*cPoise*etawat
479 if(me.eq.king.or..not.out1file)
480 & write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
482 gamp=scal_fric*(pstok+rwat)*eta
483 stdfp=dsqrt(2*Rb*t_bath/d_time)
485 gamsc(i)=scal_fric*(restok(i)+rwat)*eta
486 stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
488 if(me.eq.king.or..not.out1file)then
489 write (iout,'(/2a/)')
490 & "Radii of site types and friction coefficients and std's of",
491 & " stochastic forces of fully exposed sites"
492 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
494 write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
495 & gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
499 if(me.eq.king.or..not.out1file)then
500 write (iout,'(a)') "Berendsen bath calculation"
501 write (iout,'(a60,f10.5)') "Temperature:",t_bath
502 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
504 & write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
505 & count_reset_moment," steps"
507 & write (iout,'(a,i10,a)')
508 & "Velocities will be reset at random every",count_reset_vel,
512 if(me.eq.king.or..not.out1file)
513 & write (iout,'(a31)') "Microcanonical mode calculation"
515 if(me.eq.king.or..not.out1file)then
516 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
518 write(iout,*) "MD running with constraints."
519 write(iout,*) "Equilibration time ", eq_time, " mtus."
520 write(iout,*) "Constraining ", nfrag," fragments."
521 write(iout,*) "Length of each fragment, weight and q0:"
523 write (iout,*) "Set of restraints #",iset
525 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
526 & ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
528 write(iout,*) "constraints between ", npair, "fragments."
529 write(iout,*) "constraint pairs, weights and q0:"
531 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
532 & ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
534 write(iout,*) "angle constraints within ", nfrag_back,
535 & "backbone fragments."
536 write(iout,*) "fragment, weights:"
538 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
539 & ifrag_back(2,i,iset),wfrag_back(1,i,iset),
540 & wfrag_back(2,i,iset),wfrag_back(3,i,iset)
543 iset=mod(kolor,nset)+1
546 if(me.eq.king.or..not.out1file)
547 & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
550 c------------------------------------------------------------------------------
553 C Read molecular data.
555 implicit real*8 (a-h,o-z)
561 include 'COMMON.IOUNITS'
564 include 'COMMON.INTERACT'
565 include 'COMMON.LOCAL'
566 include 'COMMON.NAMES'
567 include 'COMMON.CHAIN'
568 include 'COMMON.FFIELD'
569 include 'COMMON.SBRIDGE'
570 include 'COMMON.HEADER'
571 include 'COMMON.CONTROL'
572 include 'COMMON.DBASE'
573 include 'COMMON.THREAD'
574 include 'COMMON.CONTACTS'
575 include 'COMMON.TORCNSTR'
576 include 'COMMON.TIME1'
577 include 'COMMON.BOUNDS'
579 include 'COMMON.SETUP'
580 character*4 sequence(maxres)
582 double precision x(maxvar)
583 character*256 pdbfile
584 character*320 weightcard
585 character*80 weightcard_t,ucase
586 dimension itype_pdb(maxres)
587 common /pizda/ itype_pdb
588 logical seq_comp,fail
589 double precision energia(0:n_ene)
595 C Read weights of the subsequent energy terms.
596 call card_concat(weightcard)
597 call reada(weightcard,'WLONG',wlong,1.0D0)
598 call reada(weightcard,'WSC',wsc,wlong)
599 call reada(weightcard,'WSCP',wscp,wlong)
600 call reada(weightcard,'WELEC',welec,1.0D0)
601 call reada(weightcard,'WVDWPP',wvdwpp,welec)
602 call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
603 call reada(weightcard,'WCORR4',wcorr4,0.0D0)
604 call reada(weightcard,'WCORR5',wcorr5,0.0D0)
605 call reada(weightcard,'WCORR6',wcorr6,0.0D0)
606 call reada(weightcard,'WTURN3',wturn3,1.0D0)
607 call reada(weightcard,'WTURN4',wturn4,1.0D0)
608 call reada(weightcard,'WTURN6',wturn6,1.0D0)
609 call reada(weightcard,'WSCCOR',wsccor,1.0D0)
610 call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
611 call reada(weightcard,'WBOND',wbond,1.0D0)
612 call reada(weightcard,'WTOR',wtor,1.0D0)
613 call reada(weightcard,'WTORD',wtor_d,1.0D0)
614 call reada(weightcard,'WANG',wang,1.0D0)
615 call reada(weightcard,'WSCLOC',wscloc,1.0D0)
616 call reada(weightcard,'SCAL14',scal14,0.4D0)
617 call reada(weightcard,'SCALSCP',scalscp,1.0d0)
618 call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
619 call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
620 call reada(weightcard,'TEMP0',temp0,300.0d0)
621 call reada(weightcard,'WLT',wliptran,0.0D0)
622 if (index(weightcard,'SOFT').gt.0) ipot=6
623 C 12/1/95 Added weight for the multi-body term WCORR
624 call reada(weightcard,'WCORRH',wcorr,1.0D0)
625 if (wcorr4.gt.0.0d0) wcorr=wcorr4
645 if(me.eq.king.or..not.out1file)
646 & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
647 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
649 10 format (/'Energy-term weights (unscaled):'//
650 & 'WSCC= ',f10.6,' (SC-SC)'/
651 & 'WSCP= ',f10.6,' (SC-p)'/
652 & 'WELEC= ',f10.6,' (p-p electr)'/
653 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
654 & 'WBOND= ',f10.6,' (stretching)'/
655 & 'WANG= ',f10.6,' (bending)'/
656 & 'WSCLOC= ',f10.6,' (SC local)'/
657 & 'WTOR= ',f10.6,' (torsional)'/
658 & 'WTORD= ',f10.6,' (double torsional)'/
659 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
660 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
661 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
662 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
663 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
664 & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
665 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
666 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
667 & 'WTURN6= ',f10.6,' (turns, 6th order)')
668 if(me.eq.king.or..not.out1file)then
669 if (wcorr4.gt.0.0d0) then
670 write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
671 & 'between contact pairs of peptide groups'
672 write (iout,'(2(a,f5.3/))')
673 & 'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
674 & 'Range of quenching the correlation terms:',2*delt_corr
675 else if (wcorr.gt.0.0d0) then
676 write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
677 & 'between contact pairs of peptide groups'
679 write (iout,'(a,f8.3)')
680 & 'Scaling factor of 1,4 SC-p interactions:',scal14
681 write (iout,'(a,f8.3)')
682 & 'General scaling factor of SC-p interactions:',scalscp
684 r0_corr=cutoff_corr-delt_corr
686 aad(i,1)=scalscp*aad(i,1)
687 aad(i,2)=scalscp*aad(i,2)
688 bad(i,1)=scalscp*bad(i,1)
689 bad(i,2)=scalscp*bad(i,2)
691 call rescale_weights(t_bath)
692 if(me.eq.king.or..not.out1file)
693 & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
694 & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
696 22 format (/'Energy-term weights (scaled):'//
697 & 'WSCC= ',f10.6,' (SC-SC)'/
698 & 'WSCP= ',f10.6,' (SC-p)'/
699 & 'WELEC= ',f10.6,' (p-p electr)'/
700 & 'WVDWPP= ',f10.6,' (p-p VDW)'/
701 & 'WBOND= ',f10.6,' (stretching)'/
702 & 'WANG= ',f10.6,' (bending)'/
703 & 'WSCLOC= ',f10.6,' (SC local)'/
704 & 'WTOR= ',f10.6,' (torsional)'/
705 & 'WTORD= ',f10.6,' (double torsional)'/
706 & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
707 & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
708 & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
709 & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
710 & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
711 & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
712 & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
713 & 'WTURN4= ',f10.6,' (turns, 4th order)'/
714 & 'WTURN6= ',f10.6,' (turns, 6th order)')
715 if(me.eq.king.or..not.out1file)
716 & write (iout,*) "Reference temperature for weights calculation:",
718 call reada(weightcard,"D0CM",d0cm,3.78d0)
719 call reada(weightcard,"AKCM",akcm,15.1d0)
720 call reada(weightcard,"AKTH",akth,11.0d0)
721 call reada(weightcard,"AKCT",akct,12.0d0)
722 call reada(weightcard,"V1SS",v1ss,-1.08d0)
723 call reada(weightcard,"V2SS",v2ss,7.61d0)
724 call reada(weightcard,"V3SS",v3ss,13.7d0)
725 call reada(weightcard,"EBR",ebr,-5.50D0)
726 dyn_ss=(index(weightcard,'DYN_SS').gt.0)
728 dyn_ss_mask(i)=.false.
732 dyn_ssbond_ij(i,j)=1.0d300
735 call reada(weightcard,"HT",Ht,0.0D0)
737 ss_depth=ebr/wsc-0.25*eps(1,1)
738 Ht=Ht/wsc-0.25*eps(1,1)
739 akcm=akcm*wstrain/wsc
740 akth=akth*wstrain/wsc
741 akct=akct*wstrain/wsc
742 v1ss=v1ss*wstrain/wsc
743 v2ss=v2ss*wstrain/wsc
744 v3ss=v3ss*wstrain/wsc
746 ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
749 if(me.eq.king.or..not.out1file) then
750 write (iout,*) "Parameters of the SS-bond potential:"
751 write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
753 write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
754 write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
755 write (iout,*)" HT",Ht
756 print *,'indpdb=',indpdb,' pdbref=',pdbref
758 if (indpdb.gt.0 .or. pdbref) then
759 read(inp,'(a)') pdbfile
760 if(me.eq.king.or..not.out1file)
761 & write (iout,'(2a)') 'PDB data will be read from file ',
762 & pdbfile(:ilen(pdbfile))
763 open(ipdbin,file=pdbfile,status='old',err=33)
765 33 write (iout,'(a)') 'Error opening PDB file.'
768 c print *,'Begin reading pdb data'
777 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3),
778 & (crefjlee(j,i+nres),j=1,3)
781 c print *,'Finished reading pdb data'
782 if(me.eq.king.or..not.out1file)
783 & write (iout,'(a,i3,a,i3)')'nsup=',nsup,
784 & ' nstart_sup=',nstart_sup
786 itype_pdb(i)=itype(i)
790 nct=nstart_sup+nsup-1
791 call contact(.false.,ncont_ref,icont_ref,co)
794 if(me.eq.king.or..not.out1file)
795 & write(iout,*)'Adding sidechains'
799 if (iti.ne.10 .and. itype(i).ne.ntyp1) then
802 do while (fail.and.nsi.le.maxsi)
803 call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
806 if(fail) write(iout,*)'Adding sidechain failed for res ',
807 & i,' after ',nsi,' trials'
810 C 10/03/12 Adam: Recalculate coordinates with new side chain positions
814 if (indpdb.eq.0) then
815 C Read sequence if not taken from the pdb file.
817 c print *,'nres=',nres
818 if (iscode.gt.0) then
819 read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
821 read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
823 C Convert sequence to numeric code
825 itype(i)=rescode(i,sequence(i),iscode)
827 C Assign initial virtual bond lengths
833 vbld(i+nres)=dsc(iabs(itype(i)))
834 vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
835 c write (iout,*) "i",i," itype",itype(i),
836 c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
840 c print '(20i4)',(itype(i),i=1,nres)
843 if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
845 if (itype(i).eq.ntyp1) then
849 else if (iabs(itype(i+1)).ne.20) then
851 else if (iabs(itype(i)).ne.20) then
858 if(me.eq.king.or..not.out1file)then
859 write (iout,*) "ITEL"
861 write (iout,*) i,itype(i),itel(i)
863 print *,'Call Read_Bridge.'
866 C 8/13/98 Set limits to generating the dihedral angles
871 read (inp,*) ndih_constr
872 if (ndih_constr.gt.0) then
874 read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
875 if(me.eq.king.or..not.out1file)then
877 & 'There are',ndih_constr,' constraints on phi angles.'
879 write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
883 phi0(i)=deg2rad*phi0(i)
884 drange(i)=deg2rad*drange(i)
886 if(me.eq.king.or..not.out1file)
887 & write (iout,*) 'FTORS',ftors
890 phibound(1,ii) = phi0(i)-drange(i)
891 phibound(2,ii) = phi0(i)+drange(i)
898 write (iout,'(a)') 'Boundaries in phi angle sampling:'
900 write (iout,'(a3,i5,2f10.1)')
901 & restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
907 cd print *,'NNT=',NNT,' NCT=',NCT
908 if (itype(1).eq.ntyp1) nnt=2
909 if (itype(nres).eq.ntyp1) nct=nct-1
911 if(me.eq.king.or..not.out1file)
912 & write (iout,'(a,i3)') 'nsup=',nsup
914 if (nsup.le.(nct-nnt+1)) then
915 do i=0,nct-nnt+1-nsup
916 if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
922 & 'Error - sequences to be superposed do not match.'
925 do i=0,nsup-(nct-nnt+1)
926 if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
928 nstart_sup=nstart_sup+i
934 & 'Error - sequences to be superposed do not match.'
937 if (nsup.eq.0) nsup=nct-nnt
938 if (nstart_sup.eq.0) nstart_sup=nnt
939 if (nstart_seq.eq.0) nstart_seq=nnt
940 if(me.eq.king.or..not.out1file)
941 & write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
942 & ' nstart_seq=',nstart_seq
944 c--- Zscore rms -------
945 if (nz_start.eq.0) nz_start=nnt
946 if (nz_end.eq.0 .and. nsup.gt.0) then
948 else if (nz_end.eq.0) then
951 if(me.eq.king.or..not.out1file)then
952 write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
953 write (iout,*) 'IZ_SC=',iz_sc
955 c----------------------
958 if (.not.pdbref) then
959 call read_angles(inp,*38)
961 38 write (iout,'(a)') 'Error reading reference structure.'
963 call MPI_Finalize(MPI_COMM_WORLD,IERROR)
964 stop 'Error reading reference structure'
966 39 call chainbuild_extconf
968 czscore call geom_to_var(nvar,coord_exp_zs(1,1))
978 call contact(.true.,ncont_ref,icont_ref,co)
980 c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
982 if (constr_dist.gt.0) call read_dist_constr
983 c write (iout,*) "After read_dist_constr nhpb",nhpb
984 if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
985 if(me.eq.king.or..not.out1file)
986 & write (iout,*) 'Contact order:',co
988 if(me.eq.king.or..not.out1file)
989 & write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
992 icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
994 if(me.eq.king.or..not.out1file)
995 & write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
996 & icont_ref(1,i),' ',
997 & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1003 if (constr_homology.gt.0) then
1004 call read_constr_homology
1005 if (indpdb.gt.0 .or. pdbref) then
1008 c(j,i)=crefjlee(j,i)
1009 cref(j,i,1)=crefjlee(j,i)
1014 write (iout,*) "Array C"
1016 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1017 & (c(j,i+nres),j=1,3)
1019 write (iout,*) "Array Cref"
1021 write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i,1),j=1,3),
1022 & (cref(j,i+nres,1),j=1,3)
1025 call int_from_cart1(.false.)
1026 call sc_loc_geom(.false.)
1028 thetaref(i)=theta(i)
1033 dc(j,i)=c(j,i+1)-c(j,i)
1034 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1039 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1040 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1048 if (nhpb.gt.0) call hpb_partition
1049 c write (iout,*) "After read_dist_constr nhpb",nhpb
1051 if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1052 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1053 & modecalc.ne.10) then
1054 C If input structure hasn't been supplied from the PDB file read or generate
1056 if (iranconf.eq.0 .and. .not. extconf) then
1057 if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1058 & write (iout,'(a)') 'Initial geometry will be read in.'
1060 read(inp,'(8f10.5)',end=36,err=36)
1061 & ((c(l,k),l=1,3),k=1,nres),
1062 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1063 write (iout,*) "Exit READ_CART"
1064 write (iout,'(8f10.5)')
1065 & ((c(l,k),l=1,3),k=1,nres),
1066 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1067 call int_from_cart1(.true.)
1068 write (iout,*) "Finish INT_TO_CART"
1071 dc(j,i)=c(j,i+1)-c(j,i)
1072 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1076 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1078 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1079 dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1085 write (iout,*) "Calling read_ang"
1086 call read_angles(inp,*36)
1087 write (iout,*) "Calling chainbuild"
1088 call chainbuild_extconf
1091 36 write (iout,'(a)') 'Error reading angle file.'
1093 call mpi_finalize( MPI_COMM_WORLD,IERR )
1095 stop 'Error reading angle file.'
1097 else if (extconf) then
1098 if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1099 & write (iout,'(a)') 'Extended chain initial geometry.'
1101 theta(i)=90d0*deg2rad
1104 phi(i)=180d0*deg2rad
1107 alph(i)=110d0*deg2rad
1110 omeg(i)=-120d0*deg2rad
1111 if (itype(i).le.0) omeg(i)=-omeg(i)
1113 c from old chainbuild
1115 C Define the origin and orientation of the coordinate system and locate the
1116 C first three CA's and SC(2).
1120 * Build the alpha-carbon chain.
1123 call locate_next_res(i)
1126 C First and last SC must coincide with the corresponding CA.
1130 dc_norm(j,nres+1)=0.0D0
1131 dc(j,nres+nres)=0.0D0
1132 dc_norm(j,nres+nres)=0.0D0
1134 c(j,nres+nres)=c(j,nres)
1137 C Define the origin and orientation of the coordinate system and locate the
1138 C first three CA's and SC(2).
1142 * Build the alpha-carbon chain.
1145 call locate_next_res(i)
1148 C First and last SC must coincide with the corresponding CA.
1152 dc_norm(j,nres+1)=0.0D0
1153 dc(j,nres+nres)=0.0D0
1154 dc_norm(j,nres+nres)=0.0D0
1156 c(j,nres+nres)=c(j,nres)
1161 if(me.eq.king.or..not.out1file)
1162 & write (iout,'(a)') 'Random-generated initial geometry.'
1166 if (me.eq.king .or. fg_rank.eq.0 .and. (
1167 & modecalc.eq.12 .or. modecalc.eq.14) ) then
1171 call gen_rand_conf(itmp,*30)
1173 30 write (iout,*) 'Failed to generate random conformation',
1174 & ', itrial=',itrial
1175 write (*,*) 'Processor:',me,
1176 & ' Failed to generate random conformation',
1186 write (iout,'(a,i3,a)') 'Processor:',me,
1187 & ' error in generating random conformation.'
1188 write (*,'(a,i3,a)') 'Processor:',me,
1189 & ' error in generating random conformation.'
1192 call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1197 & ' error in generating random conformation.'
1202 elseif (modecalc.eq.4) then
1203 read (inp,'(a)') intinname
1204 open (intin,file=intinname,status='old',err=333)
1205 if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1206 & write (iout,'(a)') 'intinname',intinname
1207 write (*,'(a)') 'Processor',myrank,' intinname',intinname
1209 333 write (iout,'(2a)') 'Error opening angle file ',intinname
1211 call MPI_Finalize(MPI_COMM_WORLD,IERR)
1213 stop 'Error opening angle file.'
1217 C Generate distance constraints, if the PDB structure is to be regularized.
1218 if (nthread.gt.0) then
1219 call read_threadbase
1222 if (me.eq.king .or. .not. out1file)
1224 if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1225 write (iout,'(/a,i3,a)')
1226 & 'The chain contains',ns,' disulfide-bridging cysteines.'
1227 write (iout,'(20i4)') (iss(i),i=1,ns)
1229 write(iout,*)"Running with dynamic disulfide-bond formation"
1231 write (iout,'(/a/)') 'Pre-formed links are:'
1237 write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1238 & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1244 if (ns.gt.0.and.dyn_ss) then
1248 forcon(i-nss)=forcon(i)
1255 dyn_ss_mask(iss(i))=.true.
1258 if (i2ndstr.gt.0) call secstrp2dihc
1259 c call geom_to_var(nvar,x)
1260 c call etotal(energia(0))
1261 c call enerprint(energia(0))
1262 c call briefout(0,etot)
1264 cd write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1265 cd write (iout,'(a)') 'Variable list:'
1266 cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1268 if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1269 & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1270 & 'Processor',myrank,': end reading molecular data.'
1274 c--------------------------------------------------------------------------
1275 logical function seq_comp(itypea,itypeb,length)
1277 integer length,itypea(length),itypeb(length)
1280 if (itypea(i).ne.itypeb(i)) then
1288 c-----------------------------------------------------------------------------
1289 subroutine read_bridge
1290 C Read information about disulfide bridges.
1291 implicit real*8 (a-h,o-z)
1292 include 'DIMENSIONS'
1296 include 'COMMON.IOUNITS'
1297 include 'COMMON.GEO'
1298 include 'COMMON.VAR'
1299 include 'COMMON.INTERACT'
1300 include 'COMMON.LOCAL'
1301 include 'COMMON.NAMES'
1302 include 'COMMON.CHAIN'
1303 include 'COMMON.FFIELD'
1304 include 'COMMON.SBRIDGE'
1305 include 'COMMON.HEADER'
1306 include 'COMMON.CONTROL'
1307 include 'COMMON.DBASE'
1308 include 'COMMON.THREAD'
1309 include 'COMMON.TIME1'
1310 include 'COMMON.SETUP'
1311 C Read bridging residues.
1312 read (inp,*) ns,(iss(i),i=1,ns)
1314 if(me.eq.king.or..not.out1file)
1315 & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1316 C Check whether the specified bridging residues are cystines.
1318 if (itype(iss(i)).ne.1) then
1319 if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1320 & 'Do you REALLY think that the residue ',
1321 & restyp(itype(iss(i))),i,
1322 & ' can form a disulfide bridge?!!!'
1323 write (*,'(2a,i3,a)')
1324 & 'Do you REALLY think that the residue ',
1325 & restyp(itype(iss(i))),i,
1326 & ' can form a disulfide bridge?!!!'
1328 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1333 C Read preformed bridges.
1335 read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1337 & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1340 C Check if the residues involved in bridges are in the specified list of
1341 C bridging residues.
1344 if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1345 & .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1346 write (iout,'(a,i3,a)') 'Disulfide pair',i,
1347 & ' contains residues present in other pairs.'
1348 write (*,'(a,i3,a)') 'Disulfide pair',i,
1349 & ' contains residues present in other pairs.'
1351 call MPI_Finalize(MPI_COMM_WORLD,ierror)
1357 if (ihpb(i).eq.iss(j)) goto 10
1359 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1362 if (jhpb(i).eq.iss(j)) goto 20
1364 write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1370 ihpb(i)=ihpb(i)+nres
1371 jhpb(i)=jhpb(i)+nres
1377 c----------------------------------------------------------------------------
1378 subroutine read_x(kanal,*)
1379 implicit real*8 (a-h,o-z)
1380 include 'DIMENSIONS'
1381 include 'COMMON.GEO'
1382 include 'COMMON.VAR'
1383 include 'COMMON.CHAIN'
1384 include 'COMMON.IOUNITS'
1385 include 'COMMON.CONTROL'
1386 include 'COMMON.LOCAL'
1387 include 'COMMON.INTERACT'
1388 c Read coordinates from input
1390 read(kanal,'(8f10.5)',end=10,err=10)
1391 & ((c(l,k),l=1,3),k=1,nres),
1392 & ((c(l,k+nres),l=1,3),k=nnt,nct)
1395 c(j,2*nres)=c(j,nres)
1397 call int_from_cart1(.false.)
1400 dc(j,i)=c(j,i+1)-c(j,i)
1401 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1405 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1407 dc(j,i+nres)=c(j,i+nres)-c(j,i)
1408 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1416 c----------------------------------------------------------------------------
1417 subroutine read_threadbase
1418 implicit real*8 (a-h,o-z)
1419 include 'DIMENSIONS'
1420 include 'COMMON.IOUNITS'
1421 include 'COMMON.GEO'
1422 include 'COMMON.VAR'
1423 include 'COMMON.INTERACT'
1424 include 'COMMON.LOCAL'
1425 include 'COMMON.NAMES'
1426 include 'COMMON.CHAIN'
1427 include 'COMMON.FFIELD'
1428 include 'COMMON.SBRIDGE'
1429 include 'COMMON.HEADER'
1430 include 'COMMON.CONTROL'
1431 include 'COMMON.DBASE'
1432 include 'COMMON.THREAD'
1433 include 'COMMON.TIME1'
1434 C Read pattern database for threading.
1435 read (icbase,*) nseq
1437 read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1438 & nres_base(2,i),nres_base(3,i)
1439 read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1441 c write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1442 c & nres_base(2,i),nres_base(3,i)
1443 c write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1447 if (weidis.eq.0.0D0) weidis=0.1D0
1456 read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1457 write (iout,'(a,i5)') 'nexcl: ',nexcl
1458 write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1461 c------------------------------------------------------------------------------
1462 subroutine setup_var
1463 implicit real*8 (a-h,o-z)
1464 include 'DIMENSIONS'
1465 include 'COMMON.IOUNITS'
1466 include 'COMMON.GEO'
1467 include 'COMMON.VAR'
1468 include 'COMMON.INTERACT'
1469 include 'COMMON.LOCAL'
1470 include 'COMMON.NAMES'
1471 include 'COMMON.CHAIN'
1472 include 'COMMON.FFIELD'
1473 include 'COMMON.SBRIDGE'
1474 include 'COMMON.HEADER'
1475 include 'COMMON.CONTROL'
1476 include 'COMMON.DBASE'
1477 include 'COMMON.THREAD'
1478 include 'COMMON.TIME1'
1479 C Set up variable list.
1485 if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1487 ialph(i,1)=nvar+nside
1491 if (indphi.gt.0) then
1493 else if (indback.gt.0) then
1498 cd write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1501 c----------------------------------------------------------------------------
1502 subroutine gen_dist_constr
1503 C Generate CA distance constraints.
1504 implicit real*8 (a-h,o-z)
1505 include 'DIMENSIONS'
1506 include 'COMMON.IOUNITS'
1507 include 'COMMON.GEO'
1508 include 'COMMON.VAR'
1509 include 'COMMON.INTERACT'
1510 include 'COMMON.LOCAL'
1511 include 'COMMON.NAMES'
1512 include 'COMMON.CHAIN'
1513 include 'COMMON.FFIELD'
1514 include 'COMMON.SBRIDGE'
1515 include 'COMMON.HEADER'
1516 include 'COMMON.CONTROL'
1517 include 'COMMON.DBASE'
1518 include 'COMMON.THREAD'
1519 include 'COMMON.TIME1'
1520 dimension itype_pdb(maxres)
1521 common /pizda/ itype_pdb
1523 cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1524 cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1525 cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1527 do i=nstart_sup,nstart_sup+nsup-1
1528 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1529 cd & ' seq_pdb', restyp(itype_pdb(i))
1530 do j=i+2,nstart_sup+nsup-1
1532 ihpb(nhpb)=i+nstart_seq-nstart_sup
1533 jhpb(nhpb)=j+nstart_seq-nstart_sup
1535 dhpb(nhpb)=dist(i,j)
1538 cd write (iout,'(a)') 'Distance constraints:'
1543 cd if (ii.gt.nres) then
1548 cd write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1549 cd & restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1550 cd & dhpb(i),forcon(i)
1554 c----------------------------------------------------------------------------
1556 implicit real*8 (a-h,o-z)
1557 include 'DIMENSIONS'
1558 include 'COMMON.MAP'
1559 include 'COMMON.IOUNITS'
1560 character*3 angid(4) /'THE','PHI','ALP','OME'/
1561 character*80 mapcard,ucase
1563 read (inp,'(a)') mapcard
1564 mapcard=ucase(mapcard)
1565 if (index(mapcard,'PHI').gt.0) then
1567 else if (index(mapcard,'THE').gt.0) then
1569 else if (index(mapcard,'ALP').gt.0) then
1571 else if (index(mapcard,'OME').gt.0) then
1574 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1575 stop 'Error - illegal variable spec in MAP card.'
1577 call readi (mapcard,'RES1',res1(imap),0)
1578 call readi (mapcard,'RES2',res2(imap),0)
1579 if (res1(imap).eq.0) then
1580 res1(imap)=res2(imap)
1581 else if (res2(imap).eq.0) then
1582 res2(imap)=res1(imap)
1584 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1586 & 'Error - illegal definition of variable group in MAP.'
1587 stop 'Error - illegal definition of variable group in MAP.'
1589 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1590 call reada(mapcard,'TO',ang_to(imap),0.0D0)
1591 call readi(mapcard,'NSTEP',nstep(imap),0)
1592 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1594 & 'Illegal boundary and/or step size specification in MAP.'
1595 stop 'Illegal boundary and/or step size specification in MAP.'
1600 c----------------------------------------------------------------------------
1602 implicit real*8 (a-h,o-z)
1603 include 'DIMENSIONS'
1604 include 'COMMON.IOUNITS'
1605 include 'COMMON.GEO'
1606 include 'COMMON.CSA'
1607 include 'COMMON.BANK'
1608 include 'COMMON.CONTROL'
1610 character*620 mcmcard
1611 call card_concat(mcmcard)
1613 call readi(mcmcard,'NCONF',nconf,50)
1614 call readi(mcmcard,'NADD',nadd,0)
1615 call readi(mcmcard,'JSTART',jstart,1)
1616 call readi(mcmcard,'JEND',jend,1)
1617 call readi(mcmcard,'NSTMAX',nstmax,500000)
1618 call readi(mcmcard,'N0',n0,1)
1619 call readi(mcmcard,'N1',n1,6)
1620 call readi(mcmcard,'N2',n2,4)
1621 call readi(mcmcard,'N3',n3,0)
1622 call readi(mcmcard,'N4',n4,0)
1623 call readi(mcmcard,'N5',n5,0)
1624 call readi(mcmcard,'N6',n6,10)
1625 call readi(mcmcard,'N7',n7,0)
1626 call readi(mcmcard,'N8',n8,0)
1627 call readi(mcmcard,'N9',n9,0)
1628 call readi(mcmcard,'N14',n14,0)
1629 call readi(mcmcard,'N15',n15,0)
1630 call readi(mcmcard,'N16',n16,0)
1631 call readi(mcmcard,'N17',n17,0)
1632 call readi(mcmcard,'N18',n18,0)
1634 vdisulf=(index(mcmcard,'DYNSS').gt.0)
1636 call readi(mcmcard,'NDIFF',ndiff,2)
1637 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1638 call readi(mcmcard,'IS1',is1,1)
1639 call readi(mcmcard,'IS2',is2,8)
1640 call readi(mcmcard,'NRAN0',nran0,4)
1641 call readi(mcmcard,'NRAN1',nran1,2)
1642 call readi(mcmcard,'IRR',irr,1)
1643 call readi(mcmcard,'NSEED',nseed,20)
1644 call readi(mcmcard,'NTOTAL',ntotal,10000)
1645 call reada(mcmcard,'CUT1',cut1,2.0d0)
1646 call reada(mcmcard,'CUT2',cut2,5.0d0)
1647 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1648 call readi(mcmcard,'ICMAX',icmax,3)
1649 call readi(mcmcard,'IRESTART',irestart,0)
1650 c!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
1653 call reada(mcmcard,'DELE',dele,20.0d0)
1654 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1655 call readi(mcmcard,'IREF',iref,0)
1656 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1657 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1658 call readi(mcmcard,'NCONF_IN',nconf_in,0)
1659 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1660 write (iout,*) "NCONF_IN",nconf_in
1663 c----------------------------------------------------------------------------
1664 cfmc subroutine mcmfread
1665 cfmc implicit real*8 (a-h,o-z)
1666 cfmc include 'DIMENSIONS'
1667 cfmc include 'COMMON.MCMF'
1668 cfmc include 'COMMON.IOUNITS'
1669 cfmc include 'COMMON.GEO'
1670 cfmc character*80 ucase
1671 cfmc character*620 mcmcard
1672 cfmc call card_concat(mcmcard)
1674 cfmc call readi(mcmcard,'MAXRANT',maxrant,1000)
1675 cfmc write(iout,*)'MAXRANT=',maxrant
1676 cfmc call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1677 cfmc write(iout,*)'MAXFAM=',maxfam
1678 cfmc call readi(mcmcard,'NNET1',nnet1,5)
1679 cfmc write(iout,*)'NNET1=',nnet1
1680 cfmc call readi(mcmcard,'NNET2',nnet2,4)
1681 cfmc write(iout,*)'NNET2=',nnet2
1682 cfmc call readi(mcmcard,'NNET3',nnet3,4)
1683 cfmc write(iout,*)'NNET3=',nnet3
1684 cfmc call readi(mcmcard,'ILASTT',ilastt,0)
1685 cfmc write(iout,*)'ILASTT=',ilastt
1686 cfmc call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1687 cfmc write(iout,*)'MAXSTR=',maxstr
1688 cfmc maxstr_f=maxstr/maxfam
1689 cfmc write(iout,*)'MAXSTR_F=',maxstr_f
1690 cfmc call readi(mcmcard,'NMCMF',nmcmf,10)
1691 cfmc write(iout,*)'NMCMF=',nmcmf
1692 cfmc call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1693 cfmc write(iout,*)'IFOCUS=',ifocus
1694 cfmc call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1695 cfmc write(iout,*)'NLOCMCMF=',nlocmcmf
1696 cfmc call readi(mcmcard,'INTPRT',intprt,1000)
1697 cfmc write(iout,*)'INTPRT=',intprt
1698 cfmc call readi(mcmcard,'IPRT',iprt,100)
1699 cfmc write(iout,*)'IPRT=',iprt
1700 cfmc call readi(mcmcard,'IMAXTR',imaxtr,100)
1701 cfmc write(iout,*)'IMAXTR=',imaxtr
1702 cfmc call readi(mcmcard,'MAXEVEN',maxeven,1000)
1703 cfmc write(iout,*)'MAXEVEN=',maxeven
1704 cfmc call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1705 cfmc write(iout,*)'MAXEVEN1=',maxeven1
1706 cfmc call readi(mcmcard,'INIMIN',inimin,200)
1707 cfmc write(iout,*)'INIMIN=',inimin
1708 cfmc call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1709 cfmc write(iout,*)'NSTEPMCMF=',nstepmcmf
1710 cfmc call readi(mcmcard,'NTHREAD',nthread,5)
1711 cfmc write(iout,*)'NTHREAD=',nthread
1712 cfmc call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1713 cfmc write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1714 cfmc call readi(mcmcard,'MAXPERT',maxpert,9)
1715 cfmc write(iout,*)'MAXPERT=',maxpert
1716 cfmc call readi(mcmcard,'IRMSD',irmsd,1)
1717 cfmc write(iout,*)'IRMSD=',irmsd
1718 cfmc call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1719 cfmc write(iout,*)'DENEMIN=',denemin
1720 cfmc call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1721 cfmc write(iout,*)'RCUT1S=',rcut1s
1722 cfmc call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1723 cfmc write(iout,*)'RCUT1E=',rcut1e
1724 cfmc call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1725 cfmc write(iout,*)'RCUT2S=',rcut2s
1726 cfmc call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1727 cfmc write(iout,*)'RCUT2E=',rcut2e
1728 cfmc call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1729 cfmc write(iout,*)'DPERT1=',d_pert1
1730 cfmc call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1731 cfmc write(iout,*)'DPERT1A=',d_pert1a
1732 cfmc call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1733 cfmc write(iout,*)'DPERT2=',d_pert2
1734 cfmc call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1735 cfmc write(iout,*)'DPERT2A=',d_pert2a
1736 cfmc call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1737 cfmc write(iout,*)'DPERT2B=',d_pert2b
1738 cfmc call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1739 cfmc write(iout,*)'DPERT2C=',d_pert2c
1740 cfmc d_pert1=deg2rad*d_pert1
1741 cfmc d_pert1a=deg2rad*d_pert1a
1742 cfmc d_pert2=deg2rad*d_pert2
1743 cfmc d_pert2a=deg2rad*d_pert2a
1744 cfmc d_pert2b=deg2rad*d_pert2b
1745 cfmc d_pert2c=deg2rad*d_pert2c
1746 cfmc call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1747 cfmc write(iout,*)'KT_MCMF1=',kt_mcmf1
1748 cfmc call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1749 cfmc write(iout,*)'KT_MCMF2=',kt_mcmf2
1750 cfmc call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1751 cfmc write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1752 cfmc call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1753 cfmc write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1754 cfmc call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1755 cfmc write(iout,*)'RCUTINI=',rcutini
1756 cfmc call reada(mcmcard,'GRAT',grat,0.5D0)
1757 cfmc write(iout,*)'GRAT=',grat
1758 cfmc call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1759 cfmc write(iout,*)'BIAS_MCMF=',bias_mcmf
1763 c----------------------------------------------------------------------------
1765 implicit real*8 (a-h,o-z)
1766 include 'DIMENSIONS'
1767 include 'COMMON.MCM'
1768 include 'COMMON.MCE'
1769 include 'COMMON.IOUNITS'
1771 character*320 mcmcard
1772 call card_concat(mcmcard)
1773 call readi(mcmcard,'MAXACC',maxacc,100)
1774 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1775 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1776 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1777 call readi(mcmcard,'MAXREPM',maxrepm,200)
1778 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1779 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1780 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1781 call reada(mcmcard,'E_UP',e_up,5.0D0)
1782 call reada(mcmcard,'DELTE',delte,0.1D0)
1783 call readi(mcmcard,'NSWEEP',nsweep,5)
1784 call readi(mcmcard,'NSTEPH',nsteph,0)
1785 call readi(mcmcard,'NSTEPC',nstepc,0)
1786 call reada(mcmcard,'TMIN',tmin,298.0D0)
1787 call reada(mcmcard,'TMAX',tmax,298.0D0)
1788 call readi(mcmcard,'NWINDOW',nwindow,0)
1789 call readi(mcmcard,'PRINT_MC',print_mc,0)
1790 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1791 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1792 ent_read=(index(mcmcard,'ENT_READ').gt.0)
1793 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1794 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1795 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1796 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1797 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1798 if (nwindow.gt.0) then
1799 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1801 winlen(i)=winend(i)-winstart(i)+1
1804 if (tmax.lt.tmin) tmax=tmin
1805 if (tmax.eq.tmin) then
1809 if (nstepc.gt.0 .and. nsteph.gt.0) then
1810 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1811 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1813 C Probabilities of different move types
1814 sumpro_type(0)=0.0D0
1815 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1816 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1817 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1818 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
1819 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1820 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1821 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1823 print *,'i',i,' sumprotype',sumpro_type(i)
1824 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1825 print *,'i',i,' sumprotype',sumpro_type(i)
1829 c----------------------------------------------------------------------------
1830 subroutine read_minim
1831 implicit real*8 (a-h,o-z)
1832 include 'DIMENSIONS'
1833 include 'COMMON.MINIM'
1834 include 'COMMON.IOUNITS'
1835 include 'COMMON.CONTROL'
1836 include 'COMMON.SETUP'
1838 character*320 minimcard
1839 call card_concat(minimcard)
1840 call readi(minimcard,'MAXMIN',maxmin,2000)
1841 call readi(minimcard,'MAXFUN',maxfun,5000)
1842 call readi(minimcard,'MINMIN',minmin,maxmin)
1843 call readi(minimcard,'MINFUN',minfun,maxmin)
1844 call reada(minimcard,'TOLF',tolf,1.0D-2)
1845 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1846 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1847 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1848 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1850 if (.not. out1file .or. me.eq.king) then
1852 write (iout,'(/80(1h*)/20x,a/80(1h*))')
1853 & 'Options in energy minimization:'
1854 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1855 & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1856 & 'MinMin:',MinMin,' MinFun:',MinFun,
1857 & ' TolF:',TolF,' RTolF:',RTolF
1863 c----------------------------------------------------------------------------
1864 subroutine read_angles(kanal,*)
1865 implicit real*8 (a-h,o-z)
1866 include 'DIMENSIONS'
1867 include 'COMMON.GEO'
1868 include 'COMMON.VAR'
1869 include 'COMMON.CHAIN'
1870 include 'COMMON.IOUNITS'
1871 include 'COMMON.CONTROL'
1872 c Read angles from input
1874 read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1875 read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1876 read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1877 read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1880 c 9/7/01 avoid 180 deg valence angle
1881 if (theta(i).gt.179.99d0) theta(i)=179.99d0
1883 theta(i)=deg2rad*theta(i)
1884 phi(i)=deg2rad*phi(i)
1885 alph(i)=deg2rad*alph(i)
1886 omeg(i)=deg2rad*omeg(i)
1891 c----------------------------------------------------------------------------
1892 subroutine reada(rekord,lancuch,wartosc,default)
1894 character*(*) rekord,lancuch
1895 double precision wartosc,default
1898 iread=index(rekord,lancuch)
1899 if (iread.eq.0) then
1903 iread=iread+ilen(lancuch)+1
1904 read (rekord(iread:),*,err=10,end=10) wartosc
1909 c----------------------------------------------------------------------------
1910 subroutine readi(rekord,lancuch,wartosc,default)
1912 character*(*) rekord,lancuch
1913 integer wartosc,default
1916 iread=index(rekord,lancuch)
1917 if (iread.eq.0) then
1921 iread=iread+ilen(lancuch)+1
1922 read (rekord(iread:),*,err=10,end=10) wartosc
1927 c----------------------------------------------------------------------------
1928 subroutine multreadi(rekord,lancuch,tablica,dim,default)
1931 integer tablica(dim),default
1932 character*(*) rekord,lancuch
1939 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1940 if (iread.eq.0) return
1941 iread=iread+ilen(lancuch)+1
1942 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1945 c----------------------------------------------------------------------------
1946 subroutine multreada(rekord,lancuch,tablica,dim,default)
1949 double precision tablica(dim),default
1950 character*(*) rekord,lancuch
1957 iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1958 if (iread.eq.0) return
1959 iread=iread+ilen(lancuch)+1
1960 read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1963 c----------------------------------------------------------------------------
1964 subroutine openunits
1965 implicit real*8 (a-h,o-z)
1966 include 'DIMENSIONS'
1969 character*16 form,nodename
1972 include 'COMMON.SETUP'
1973 include 'COMMON.IOUNITS'
1975 include 'COMMON.CONTROL'
1976 integer lenpre,lenpot,ilen,lentmp
1978 character*3 out1file_text,ucase
1981 c print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1982 call getenv_loc("PREFIX",prefix)
1984 call getenv_loc("POT",pot)
1985 call getenv_loc("DIRTMP",tmpdir)
1986 call getenv_loc("CURDIR",curdir)
1987 call getenv_loc("OUT1FILE",out1file_text)
1988 c print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1989 out1file_text=ucase(out1file_text)
1990 if (out1file_text(1:1).eq."Y") then
1993 out1file=fg_rank.gt.0
1998 if (lentmp.gt.0) then
1999 write (*,'(80(1h!))')
2000 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
2001 write (*,'(80(1h!))')
2002 write (*,*)"All output files will be on node /tmp directory."
2004 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2005 if (me.eq.king) then
2006 write (*,*) "The master node is ",nodename
2007 else if (fg_rank.eq.0) then
2008 write (*,*) "I am the CG slave node ",nodename
2010 write (*,*) "I am the FG slave node ",nodename
2013 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2014 lenpre = lentmp+lenpre+1
2016 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2017 C Get the names and open the input files
2018 #if defined(WINIFL) || defined(WINPGI)
2019 open(1,file=pref_orig(:ilen(pref_orig))//
2020 & '.inp',status='old',readonly,shared)
2021 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2022 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2023 C Get parameter filenames and open the parameter files.
2024 call getenv_loc('BONDPAR',bondname)
2025 open (ibond,file=bondname,status='old',readonly,shared)
2026 call getenv_loc('THETPAR',thetname)
2027 open (ithep,file=thetname,status='old',readonly,shared)
2029 call getenv_loc('THETPARPDB',thetname_pdb)
2030 open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2032 call getenv_loc('ROTPAR',rotname)
2033 open (irotam,file=rotname,status='old',readonly,shared)
2035 call getenv_loc('ROTPARPDB',rotname_pdb)
2036 open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2038 call getenv_loc('TORPAR',torname)
2039 open (itorp,file=torname,status='old',readonly,shared)
2040 call getenv_loc('TORDPAR',tordname)
2041 open (itordp,file=tordname,status='old',readonly,shared)
2042 call getenv_loc('FOURIER',fouriername)
2043 open (ifourier,file=fouriername,status='old',readonly,shared)
2044 call getenv_loc('ELEPAR',elename)
2045 open (ielep,file=elename,status='old',readonly,shared)
2046 call getenv_loc('SIDEPAR',sidename)
2047 open (isidep,file=sidename,status='old',readonly,shared)
2048 #elif (defined CRAY) || (defined AIX)
2049 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2051 c print *,"Processor",myrank," opened file 1"
2052 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2053 c print *,"Processor",myrank," opened file 9"
2054 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2055 C Get parameter filenames and open the parameter files.
2056 call getenv_loc('BONDPAR',bondname)
2057 open (ibond,file=bondname,status='old',action='read')
2058 c print *,"Processor",myrank," opened file IBOND"
2059 call getenv_loc('THETPAR',thetname)
2060 open (ithep,file=thetname,status='old',action='read')
2061 c print *,"Processor",myrank," opened file ITHEP"
2063 call getenv_loc('THETPARPDB',thetname_pdb)
2064 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2066 call getenv_loc('ROTPAR',rotname)
2067 open (irotam,file=rotname,status='old',action='read')
2068 c print *,"Processor",myrank," opened file IROTAM"
2070 call getenv_loc('ROTPARPDB',rotname_pdb)
2071 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2073 call getenv_loc('TORPAR',torname)
2074 open (itorp,file=torname,status='old',action='read')
2075 c print *,"Processor",myrank," opened file ITORP"
2076 call getenv_loc('TORDPAR',tordname)
2077 open (itordp,file=tordname,status='old',action='read')
2078 c print *,"Processor",myrank," opened file ITORDP"
2079 call getenv_loc('SCCORPAR',sccorname)
2080 open (isccor,file=sccorname,status='old',action='read')
2081 c print *,"Processor",myrank," opened file ISCCOR"
2082 call getenv_loc('FOURIER',fouriername)
2083 open (ifourier,file=fouriername,status='old',action='read')
2084 c print *,"Processor",myrank," opened file IFOURIER"
2085 call getenv_loc('ELEPAR',elename)
2086 open (ielep,file=elename,status='old',action='read')
2087 c print *,"Processor",myrank," opened file IELEP"
2088 call getenv_loc('SIDEPAR',sidename)
2089 open (isidep,file=sidename,status='old',action='read')
2090 c print *,"Processor",myrank," opened file ISIDEP"
2091 c print *,"Processor",myrank," opened parameter files"
2093 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2094 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2095 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2096 C Get parameter filenames and open the parameter files.
2097 call getenv_loc('BONDPAR',bondname)
2098 open (ibond,file=bondname,status='old')
2099 call getenv_loc('THETPAR',thetname)
2100 open (ithep,file=thetname,status='old')
2102 call getenv_loc('THETPARPDB',thetname_pdb)
2103 open (ithep_pdb,file=thetname_pdb,status='old')
2105 call getenv_loc('ROTPAR',rotname)
2106 open (irotam,file=rotname,status='old')
2108 call getenv_loc('ROTPARPDB',rotname_pdb)
2109 open (irotam_pdb,file=rotname_pdb,status='old')
2111 call getenv_loc('TORPAR',torname)
2112 open (itorp,file=torname,status='old')
2113 call getenv_loc('TORDPAR',tordname)
2114 open (itordp,file=tordname,status='old')
2115 call getenv_loc('SCCORPAR',sccorname)
2116 open (isccor,file=sccorname,status='old')
2117 call getenv_loc('FOURIER',fouriername)
2118 open (ifourier,file=fouriername,status='old')
2119 call getenv_loc('ELEPAR',elename)
2120 open (ielep,file=elename,status='old')
2121 call getenv_loc('SIDEPAR',sidename)
2122 open (isidep,file=sidename,status='old')
2124 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2126 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2127 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2128 C Get parameter filenames and open the parameter files.
2129 call getenv_loc('BONDPAR',bondname)
2130 open (ibond,file=bondname,status='old',readonly)
2131 call getenv_loc('THETPAR',thetname)
2132 open (ithep,file=thetname,status='old',readonly)
2133 call getenv_loc('ROTPAR',rotname)
2134 open (irotam,file=rotname,status='old',readonly)
2135 call getenv_loc('TORPAR',torname)
2136 open (itorp,file=torname,status='old',readonly)
2137 call getenv_loc('TORDPAR',tordname)
2138 open (itordp,file=tordname,status='old',readonly)
2139 call getenv_loc('SCCORPAR',sccorname)
2140 open (isccor,file=sccorname,status='old',readonly)
2142 call getenv_loc('THETPARPDB',thetname_pdb)
2143 c print *,"thetname_pdb ",thetname_pdb
2144 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2145 c print *,ithep_pdb," opened"
2147 call getenv_loc('FOURIER',fouriername)
2148 open (ifourier,file=fouriername,status='old',readonly)
2149 call getenv_loc('ELEPAR',elename)
2150 open (ielep,file=elename,status='old',readonly)
2151 call getenv_loc('SIDEPAR',sidename)
2152 open (isidep,file=sidename,status='old',readonly)
2153 call getenv_loc('LIPTRANPAR',liptranname)
2154 open (iliptranpar,file=liptranname,status='old',action='read')
2156 call getenv_loc('ROTPARPDB',rotname_pdb)
2157 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2162 C 8/9/01 In the newest version SCp interaction constants are read from a file
2163 C Use -DOLDSCP to use hard-coded constants instead.
2165 call getenv_loc('SCPPAR',scpname)
2166 #if defined(WINIFL) || defined(WINPGI)
2167 open (iscpp,file=scpname,status='old',readonly,shared)
2168 #elif (defined CRAY) || (defined AIX)
2169 open (iscpp,file=scpname,status='old',action='read')
2171 open (iscpp,file=scpname,status='old')
2173 open (iscpp,file=scpname,status='old',readonly)
2176 call getenv_loc('PATTERN',patname)
2177 #if defined(WINIFL) || defined(WINPGI)
2178 open (icbase,file=patname,status='old',readonly,shared)
2179 #elif (defined CRAY) || (defined AIX)
2180 open (icbase,file=patname,status='old',action='read')
2182 open (icbase,file=patname,status='old')
2184 open (icbase,file=patname,status='old',readonly)
2187 C Open output file only for CG processes
2188 c print *,"Processor",myrank," fg_rank",fg_rank
2189 if (fg_rank.eq.0) then
2191 if (nodes.eq.1) then
2194 npos = dlog10(dfloat(nodes-1))+1
2196 if (npos.lt.3) npos=3
2197 write (liczba,'(i1)') npos
2198 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2200 write (liczba,form) me
2201 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2202 & liczba(:ilen(liczba))
2203 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2205 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2207 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2208 & liczba(:ilen(liczba))//'.mol2'
2209 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2210 & liczba(:ilen(liczba))//'.stat'
2212 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2213 & //liczba(:ilen(liczba))//'.stat')
2214 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2217 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2218 & liczba(:ilen(liczba))//'.const'
2223 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2224 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2225 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2226 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2227 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2229 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2231 rest2name=prefix(:ilen(prefix))//'.rst'
2233 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2236 #if defined(AIX) || defined(PGI)
2237 if (me.eq.king .or. .not. out1file)
2238 & open(iout,file=outname,status='unknown')
2240 if (fg_rank.gt.0) then
2241 write (liczba,'(i3.3)') myrank/nfgtasks
2242 write (ll,'(bz,i3.3)') fg_rank
2243 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2248 open(igeom,file=intname,status='unknown',position='append')
2249 open(ipdb,file=pdbname,status='unknown')
2250 open(imol2,file=mol2name,status='unknown')
2251 open(istat,file=statname,status='unknown',position='append')
2253 c1out open(iout,file=outname,status='unknown')
2256 if (me.eq.king .or. .not.out1file)
2257 & open(iout,file=outname,status='unknown')
2259 if (fg_rank.gt.0) then
2260 write (liczba,'(i3.3)') myrank/nfgtasks
2261 write (ll,'(bz,i3.3)') fg_rank
2262 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2267 open(igeom,file=intname,status='unknown',access='append')
2268 open(ipdb,file=pdbname,status='unknown')
2269 open(imol2,file=mol2name,status='unknown')
2270 open(istat,file=statname,status='unknown',access='append')
2272 c1out open(iout,file=outname,status='unknown')
2275 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2276 csa_seed=prefix(:lenpre)//'.CSA.seed'
2277 csa_history=prefix(:lenpre)//'.CSA.history'
2278 csa_bank=prefix(:lenpre)//'.CSA.bank'
2279 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2280 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2281 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2282 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2283 csa_int=prefix(:lenpre)//'.int'
2284 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2285 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2286 csa_in=prefix(:lenpre)//'.CSA.in'
2287 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2290 write (iout,'(80(1h-))')
2291 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2292 write (iout,'(80(1h-))')
2293 write (iout,*) "Input file : ",
2294 & pref_orig(:ilen(pref_orig))//'.inp'
2295 write (iout,*) "Output file : ",
2296 & outname(:ilen(outname))
2298 write (iout,*) "Sidechain potential file : ",
2299 & sidename(:ilen(sidename))
2301 write (iout,*) "SCp potential file : ",
2302 & scpname(:ilen(scpname))
2304 write (iout,*) "Electrostatic potential file : ",
2305 & elename(:ilen(elename))
2306 write (iout,*) "Cumulant coefficient file : ",
2307 & fouriername(:ilen(fouriername))
2308 write (iout,*) "Torsional parameter file : ",
2309 & torname(:ilen(torname))
2310 write (iout,*) "Double torsional parameter file : ",
2311 & tordname(:ilen(tordname))
2312 write (iout,*) "SCCOR parameter file : ",
2313 & sccorname(:ilen(sccorname))
2314 write (iout,*) "Bond & inertia constant file : ",
2315 & bondname(:ilen(bondname))
2316 write (iout,*) "Bending parameter file : ",
2317 & thetname(:ilen(thetname))
2318 write (iout,*) "Rotamer parameter file : ",
2319 & rotname(:ilen(rotname))
2320 write (iout,*) "Thetpdb parameter file : ",
2321 & thetname_pdb(:ilen(thetname_pdb))
2322 write (iout,*) "Threading database : ",
2323 & patname(:ilen(patname))
2325 &write (iout,*)" DIRTMP : ",
2327 write (iout,'(80(1h-))')
2331 c----------------------------------------------------------------------------
2332 subroutine card_concat(card)
2333 implicit real*8 (a-h,o-z)
2334 include 'DIMENSIONS'
2335 include 'COMMON.IOUNITS'
2337 character*80 karta,ucase
2339 read (inp,'(a)') karta
2342 do while (karta(80:80).eq.'&')
2343 card=card(:ilen(card)+1)//karta(:79)
2344 read (inp,'(a)') karta
2347 card=card(:ilen(card)+1)//karta
2350 c----------------------------------------------------------------------------------
2352 implicit real*8 (a-h,o-z)
2353 include 'DIMENSIONS'
2354 include 'COMMON.CHAIN'
2355 include 'COMMON.IOUNITS'
2357 include 'COMMON.CONTROL'
2358 open(irest2,file=rest2name,status='unknown')
2359 read(irest2,*) totT,EK,potE,totE,t_bath
2362 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2365 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2367 if(usampl.or.homol_nset.gt.1) then
2368 read (irest2,*) iset
2373 c---------------------------------------------------------------------------------
2374 subroutine read_fragments
2375 implicit real*8 (a-h,o-z)
2376 include 'DIMENSIONS'
2380 include 'COMMON.SETUP'
2381 include 'COMMON.CHAIN'
2382 include 'COMMON.IOUNITS'
2384 include 'COMMON.CONTROL'
2386 read(inp,*) nset,nfrag,npair,nfrag_back
2387 if(me.eq.king.or..not.out1file)
2388 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2389 & " nfrag_back",nfrag_back
2391 read(inp,*) mset(iset1)
2393 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2395 if(me.eq.king.or..not.out1file)
2396 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2397 & ifrag(2,i,iset1), qinfrag(i,iset1)
2400 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2402 if(me.eq.king.or..not.out1file)
2403 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2404 & ipair(2,i,iset1), qinpair(i,iset1)
2407 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2408 & wfrag_back(3,i,iset1),
2409 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2410 if(me.eq.king.or..not.out1file)
2411 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2412 & wfrag_back(2,i,iset1),
2413 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2414 & ifrag_back(2,i,iset1)
2419 C---------------------------------------------------------------------------
2420 subroutine read_afminp
2421 implicit real*8 (a-h,o-z)
2422 include 'DIMENSIONS'
2426 include 'COMMON.SETUP'
2427 include 'COMMON.CONTROL'
2428 include 'COMMON.CHAIN'
2429 include 'COMMON.IOUNITS'
2430 include 'COMMON.SBRIDGE'
2431 character*320 afmcard
2432 c print *, "wchodze"
2433 call card_concat(afmcard)
2434 call readi(afmcard,"BEG",afmbeg,0)
2435 call readi(afmcard,"END",afmend,0)
2436 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2437 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2438 print *,'FORCE=' ,forceAFMconst
2439 CCCC NOW PROPERTIES FOR AFM
2442 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2444 distafminit=dsqrt(distafminit)
2445 print *,'initdist',distafminit
2449 c-------------------------------------------------------------------------------
2450 subroutine read_dist_constr
2451 implicit real*8 (a-h,o-z)
2452 include 'DIMENSIONS'
2456 include 'COMMON.SETUP'
2457 include 'COMMON.CONTROL'
2458 include 'COMMON.CHAIN'
2459 include 'COMMON.IOUNITS'
2460 include 'COMMON.SBRIDGE'
2461 integer ifrag_(2,100),ipair_(2,100)
2462 double precision wfrag_(100),wpair_(100)
2463 character*500 controlcard
2464 c write (iout,*) "Calling read_dist_constr"
2465 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2467 call card_concat(controlcard)
2468 call readi(controlcard,"NFRAG",nfrag_,0)
2469 call readi(controlcard,"NPAIR",npair_,0)
2470 call readi(controlcard,"NDIST",ndist_,0)
2471 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2472 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2473 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2474 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2475 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2476 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2477 c write (iout,*) "IFRAG"
2479 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2481 c write (iout,*) "IPAIR"
2483 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2485 if (.not.refstr .and. nfrag.gt.0) then
2487 & "ERROR: no reference structure to compute distance restraints"
2489 & "Restraints must be specified explicitly (NDIST=number)"
2492 if (nfrag.lt.2 .and. npair.gt.0) then
2493 write (iout,*) "ERROR: Less than 2 fragments specified",
2494 & " but distance restraints between pairs requested"
2499 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2500 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2501 & ifrag_(2,i)=nstart_sup+nsup-1
2502 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2504 if (wfrag_(i).gt.0.0d0) then
2505 do j=ifrag_(1,i),ifrag_(2,i)-1
2506 do k=j+1,ifrag_(2,i)
2507 c write (iout,*) "j",j," k",k
2509 if (constr_dist.eq.1) then
2514 forcon(nhpb)=wfrag_(i)
2515 else if (constr_dist.eq.2) then
2516 if (ddjk.le.dist_cut) then
2521 forcon(nhpb)=wfrag_(i)
2528 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2531 if (.not.out1file .or. me.eq.king)
2532 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2533 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2535 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2536 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2543 if (wpair_(i).gt.0.0d0) then
2551 do j=ifrag_(1,ii),ifrag_(2,ii)
2552 do k=ifrag_(1,jj),ifrag_(2,jj)
2556 forcon(nhpb)=wpair_(i)
2557 dhpb(nhpb)=dist(j,k)
2559 if (.not.out1file .or. me.eq.king)
2560 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2561 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2563 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2564 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2571 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2572 if (forcon(nhpb+1).gt.0.0d0) then
2574 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2576 if (.not.out1file .or. me.eq.king)
2577 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2578 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2580 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2581 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2588 c-------------------------------------------------------------------------------
2590 subroutine read_constr_homology
2592 include 'DIMENSIONS'
2596 include 'COMMON.SETUP'
2597 include 'COMMON.CONTROL'
2598 include 'COMMON.CHAIN'
2599 include 'COMMON.IOUNITS'
2601 include 'COMMON.GEO'
2602 include 'COMMON.INTERACT'
2603 include 'COMMON.NAMES'
2605 c For new homol impl
2607 include 'COMMON.VAR'
2610 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2612 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2613 c & sigma_odl_temp(maxres,maxres,max_template)
2615 character*24 model_ki_dist, model_ki_angle
2616 character*500 controlcard
2617 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2618 logical lprn /.true./
2623 c FP - Nov. 2014 Temporary specifications for new vars
2625 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2627 double precision, dimension (max_template,maxres) :: rescore
2628 double precision, dimension (max_template,maxres) :: rescore2
2629 double precision, dimension (max_template,maxres) :: rescore3
2630 character*24 pdbfile,tpl_k_rescore
2631 c -----------------------------------------------------------------
2632 c Reading multiple PDB ref structures and calculation of retraints
2633 c not using pre-computed ones stored in files model_ki_{dist,angle}
2635 c -----------------------------------------------------------------
2638 c Alternative: reading from input
2639 call card_concat(controlcard)
2640 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2641 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2642 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2643 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2644 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2645 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2646 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2647 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2648 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2649 if(.not.read2sigma.and.start_from_model) then
2650 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2651 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2652 start_from_model=.false.
2654 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2655 & write(iout,*) 'START_FROM_MODELS is ON'
2656 if(start_from_model .and. rest) then
2657 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2658 write(iout,*) 'START_FROM_MODELS is OFF'
2659 write(iout,*) 'remove restart keyword from input'
2662 if (homol_nset.gt.1)then
2663 call card_concat(controlcard)
2664 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2665 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2666 write(iout,*) "iset homology_weight "
2668 write(iout,*) i,waga_homology(i)
2671 iset=mod(kolor,homol_nset)+1
2674 waga_homology(1)=1.0
2677 cd write (iout,*) "nnt",nnt," nct",nct
2684 c write(iout,*) 'nnt=',nnt,'nct=',nct
2687 do k=1,constr_homology
2700 do k=1,constr_homology
2702 read(inp,'(a)') pdbfile
2703 if(me.eq.king .or. .not. out1file)
2704 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2705 & pdbfile(:ilen(pdbfile))
2706 open(ipdbin,file=pdbfile,status='old',err=33)
2708 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2709 & pdbfile(:ilen(pdbfile))
2712 c print *,'Begin reading pdb data'
2714 c Files containing res sim or local scores (former containing sigmas)
2717 write(kic2,'(bz,i2.2)') k
2719 tpl_k_rescore="template"//kic2//".sco"
2722 if (read2sigma) then
2723 call readpdb_template(k)
2728 c Distance restraints
2731 C Copy the coordinates from reference coordinates (?)
2735 c write (iout,*) "c(",j,i,") =",c(j,i)
2739 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2741 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2742 open (ientin,file=tpl_k_rescore,status='old')
2743 if (nnt.gt.1) rescore(k,1)=0.0d0
2744 do irec=nnt,nct ! loop for reading res sim
2745 if (read2sigma) then
2746 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2747 & rescore3_tmp,idomain_tmp
2749 idomain(k,i_tmp)=idomain_tmp
2750 rescore(k,i_tmp)=rescore_tmp
2751 rescore2(k,i_tmp)=rescore2_tmp
2752 rescore3(k,i_tmp)=rescore3_tmp
2753 if (.not. out1file .or. me.eq.king)
2754 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
2755 & i_tmp,rescore2_tmp,rescore_tmp,
2756 & rescore3_tmp,idomain_tmp
2759 read (ientin,*,end=1401) rescore_tmp
2761 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2762 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2763 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2768 if (waga_dist.ne.0.0d0) then
2776 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2777 c write (iout,*) k,i,j,distal,dist2_cut
2779 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2780 & .and. distal.le.dist2_cut ) then
2786 c write (iout,*) "k",k
2787 c write (iout,*) "i",i," j",j," constr_homology",
2792 if (read2sigma) then
2795 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2797 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2798 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2799 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2801 if (odl(k,ii).le.dist_cut) then
2802 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2805 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2806 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2808 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2809 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2813 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2816 l_homo(k,ii)=.false.
2823 c Theta, dihedral and SC retraints
2825 if (waga_angle.gt.0.0d0) then
2826 c open (ientin,file=tpl_k_sigma_dih,status='old')
2827 c do irec=1,maxres-3 ! loop for reading sigma_dih
2828 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2829 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2830 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2831 c & sigma_dih(k,i+nnt-1)
2836 if (idomain(k,i).eq.0) then
2840 dih(k,i)=phiref(i) ! right?
2841 c read (ientin,*) sigma_dih(k,i) ! original variant
2842 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2843 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2844 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2845 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2846 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2848 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2849 & rescore(k,i-2)+rescore(k,i-3))/4.0
2850 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2851 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2852 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2853 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2854 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2855 c Instead of res sim other local measure of b/b str reliability possible
2856 if (sigma_dih(k,i).ne.0)
2857 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2858 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2863 if (waga_theta.gt.0.0d0) then
2864 c open (ientin,file=tpl_k_sigma_theta,status='old')
2865 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2866 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2867 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2868 c & sigma_theta(k,i+nnt-1)
2873 do i = nnt+2,nct ! right? without parallel.
2874 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2875 c do i=ithet_start,ithet_end ! with FG parallel.
2876 if (idomain(k,i).eq.0) then
2877 sigma_theta(k,i)=0.0
2880 thetatpl(k,i)=thetaref(i)
2881 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2882 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2883 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2884 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2885 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2886 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2887 & rescore(k,i-2))/3.0
2888 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2889 if (sigma_theta(k,i).ne.0)
2890 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2892 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2893 c rescore(k,i-2) ! right expression ?
2894 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2898 if (waga_d.gt.0.0d0) then
2899 c open (ientin,file=tpl_k_sigma_d,status='old')
2900 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2901 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2902 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2903 c & sigma_d(k,i+nnt-1)
2907 do i = nnt,nct ! right? without parallel.
2908 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2909 c do i=loc_start,loc_end ! with FG parallel.
2910 if (itype(i).eq.10) cycle
2911 if (idomain(k,i).eq.0 ) then
2918 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2919 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2920 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2921 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
2922 sigma_d(k,i)=rescore3(k,i) ! right expression ?
2923 if (sigma_d(k,i).ne.0)
2924 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2926 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
2927 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2928 c read (ientin,*) sigma_d(k,i) ! 1st variant
2933 c remove distance restraints not used in any model from the list
2934 c shift data in all arrays
2936 if (waga_dist.ne.0.0d0) then
2942 if (ii_in_use(ii).eq.0.and.liiflag) then
2946 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
2947 & .not.liiflag.and.ii.eq.lim_odl) then
2948 if (ii.eq.lim_odl) then
2949 iishift=ii-iistart+1
2954 do ki=iistart,lim_odl-iishift
2955 ires_homo(ki)=ires_homo(ki+iishift)
2956 jres_homo(ki)=jres_homo(ki+iishift)
2957 ii_in_use(ki)=ii_in_use(ki+iishift)
2958 do k=1,constr_homology
2959 odl(k,ki)=odl(k,ki+iishift)
2960 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
2961 l_homo(k,ki)=l_homo(k,ki+iishift)
2965 lim_odl=lim_odl-iishift
2970 if (constr_homology.gt.0) call homology_partition
2971 if (constr_homology.gt.0) call init_int_table
2972 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2973 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2977 if (.not.lprn) return
2978 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2979 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2980 write (iout,*) "Distance restraints from templates"
2982 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
2983 & ii,ires_homo(ii),jres_homo(ii),
2984 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
2985 & ki=1,constr_homology)
2987 write (iout,*) "Dihedral angle restraints from templates"
2989 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
2990 & (rad2deg*dih(ki,i),
2991 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2993 write (iout,*) "Virtual-bond angle restraints from templates"
2995 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
2996 & (rad2deg*thetatpl(ki,i),
2997 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2999 write (iout,*) "SC restraints from templates"
3001 write(iout,'(i5,100(4f8.2,4x))') i,
3002 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3003 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3006 c -----------------------------------------------------------------
3009 c----------------------------------------------------------------------
3012 subroutine flush(iu)
3017 subroutine flush(iu)
3022 c------------------------------------------------------------------------------
3023 subroutine copy_to_tmp(source)
3024 include "DIMENSIONS"
3025 include "COMMON.IOUNITS"
3026 character*(*) source
3027 character* 256 tmpfile
3031 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3032 inquire(file=tmpfile,exist=ex)
3034 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3035 & " to temporary directory..."
3036 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3037 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3041 c------------------------------------------------------------------------------
3042 subroutine move_from_tmp(source)
3043 include "DIMENSIONS"
3044 include "COMMON.IOUNITS"
3045 character*(*) source
3048 write (*,*) "Moving ",source(:ilen(source)),
3049 & " from temporary directory to working directory"
3050 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3051 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3054 c------------------------------------------------------------------------------
3055 subroutine random_init(seed)
3057 C Initialize random number generator
3059 implicit real*8 (a-h,o-z)
3060 include 'DIMENSIONS'
3066 logical OKRandom, prng_restart
3068 integer iseed_array(4)
3070 include 'COMMON.IOUNITS'
3071 include 'COMMON.TIME1'
3072 include 'COMMON.THREAD'
3073 include 'COMMON.SBRIDGE'
3074 include 'COMMON.CONTROL'
3075 include 'COMMON.MCM'
3076 include 'COMMON.MAP'
3077 include 'COMMON.HEADER'
3078 include 'COMMON.CSA'
3079 include 'COMMON.CHAIN'
3080 include 'COMMON.MUCA'
3082 include 'COMMON.FFIELD'
3083 include 'COMMON.SETUP'
3084 iseed=-dint(dabs(seed))
3085 if (iseed.eq.0) then
3086 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3087 & 'Random seed undefined. The program will stop.'
3088 write (*,'(/80(1h*)/20x,a/80(1h*))')
3089 & 'Random seed undefined. The program will stop.'
3091 call mpi_finalize(mpi_comm_world,ierr)
3093 stop 'Bad random seed.'
3096 if (fg_rank.eq.0) then
3100 if(me.eq.king .or. .not. out1file)
3101 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3102 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3103 OKRandom = prng_restart(me,iseedi8)
3106 tmp=65536.0d0**(4-i)
3107 iseed_array(i) = dint(seed/tmp)
3108 seed=seed-iseed_array(i)*tmp
3111 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3112 & (iseed_array(i),i=1,4)
3113 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3114 & (iseed_array(i),i=1,4)
3115 OKRandom = prng_restart(me,iseed_array)
3118 c r1 = prng_next(me)
3119 r1=ran_number(0.0D0,1.0D0)
3120 if(me.eq.king .or. .not. out1file)
3121 & write (iout,*) 'ran_num',r1
3122 if (r1.lt.0.0d0) OKRandom=.false.
3124 if (.not.OKRandom) then
3125 write (iout,*) 'PRNG IS NOT WORKING!!!'
3126 print *,'PRNG IS NOT WORKING!!!'
3129 call mpi_abort(mpi_comm_world,error_msg,ierr)
3132 write (iout,*) 'too many processors for parallel prng'
3133 write (*,*) 'too many processors for parallel prng'
3141 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)