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')
2123 call getenv_loc('LIPTRANPAR',liptranname)
2124 open (iliptranpar,file=liptranname,status='old')
2126 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2128 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2129 C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2130 C Get parameter filenames and open the parameter files.
2131 call getenv_loc('BONDPAR',bondname)
2132 open (ibond,file=bondname,status='old',readonly)
2133 call getenv_loc('THETPAR',thetname)
2134 open (ithep,file=thetname,status='old',readonly)
2135 call getenv_loc('ROTPAR',rotname)
2136 open (irotam,file=rotname,status='old',readonly)
2137 call getenv_loc('TORPAR',torname)
2138 open (itorp,file=torname,status='old',readonly)
2139 call getenv_loc('TORDPAR',tordname)
2140 open (itordp,file=tordname,status='old',readonly)
2141 call getenv_loc('SCCORPAR',sccorname)
2142 open (isccor,file=sccorname,status='old',readonly)
2144 call getenv_loc('THETPARPDB',thetname_pdb)
2145 c print *,"thetname_pdb ",thetname_pdb
2146 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2147 c print *,ithep_pdb," opened"
2149 call getenv_loc('FOURIER',fouriername)
2150 open (ifourier,file=fouriername,status='old',readonly)
2151 call getenv_loc('ELEPAR',elename)
2152 open (ielep,file=elename,status='old',readonly)
2153 call getenv_loc('SIDEPAR',sidename)
2154 open (isidep,file=sidename,status='old',readonly)
2155 call getenv_loc('LIPTRANPAR',liptranname)
2156 open (iliptranpar,file=liptranname,status='old',readonly)
2158 call getenv_loc('ROTPARPDB',rotname_pdb)
2159 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2164 C 8/9/01 In the newest version SCp interaction constants are read from a file
2165 C Use -DOLDSCP to use hard-coded constants instead.
2167 call getenv_loc('SCPPAR',scpname)
2168 #if defined(WINIFL) || defined(WINPGI)
2169 open (iscpp,file=scpname,status='old',readonly,shared)
2170 #elif (defined CRAY) || (defined AIX)
2171 open (iscpp,file=scpname,status='old',action='read')
2173 open (iscpp,file=scpname,status='old')
2175 open (iscpp,file=scpname,status='old',readonly)
2178 call getenv_loc('PATTERN',patname)
2179 #if defined(WINIFL) || defined(WINPGI)
2180 open (icbase,file=patname,status='old',readonly,shared)
2181 #elif (defined CRAY) || (defined AIX)
2182 open (icbase,file=patname,status='old',action='read')
2184 open (icbase,file=patname,status='old')
2186 open (icbase,file=patname,status='old',readonly)
2189 C Open output file only for CG processes
2190 c print *,"Processor",myrank," fg_rank",fg_rank
2191 if (fg_rank.eq.0) then
2193 if (nodes.eq.1) then
2196 npos = dlog10(dfloat(nodes-1))+1
2198 if (npos.lt.3) npos=3
2199 write (liczba,'(i1)') npos
2200 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2202 write (liczba,form) me
2203 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2204 & liczba(:ilen(liczba))
2205 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2207 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2209 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2210 & liczba(:ilen(liczba))//'.mol2'
2211 statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2212 & liczba(:ilen(liczba))//'.stat'
2214 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2215 & //liczba(:ilen(liczba))//'.stat')
2216 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2219 qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2220 & liczba(:ilen(liczba))//'.const'
2225 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2226 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2227 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2228 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2229 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2231 & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2233 rest2name=prefix(:ilen(prefix))//'.rst'
2235 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2238 #if defined(AIX) || defined(PGI)
2239 if (me.eq.king .or. .not. out1file)
2240 & open(iout,file=outname,status='unknown')
2242 if (fg_rank.gt.0) then
2243 write (liczba,'(i3.3)') myrank/nfgtasks
2244 write (ll,'(bz,i3.3)') fg_rank
2245 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2250 open(igeom,file=intname,status='unknown',position='append')
2251 open(ipdb,file=pdbname,status='unknown')
2252 open(imol2,file=mol2name,status='unknown')
2253 open(istat,file=statname,status='unknown',position='append')
2255 c1out open(iout,file=outname,status='unknown')
2258 if (me.eq.king .or. .not.out1file)
2259 & open(iout,file=outname,status='unknown')
2261 if (fg_rank.gt.0) then
2262 write (liczba,'(i3.3)') myrank/nfgtasks
2263 write (ll,'(bz,i3.3)') fg_rank
2264 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2269 open(igeom,file=intname,status='unknown',access='append')
2270 open(ipdb,file=pdbname,status='unknown')
2271 open(imol2,file=mol2name,status='unknown')
2272 open(istat,file=statname,status='unknown',access='append')
2274 c1out open(iout,file=outname,status='unknown')
2277 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2278 csa_seed=prefix(:lenpre)//'.CSA.seed'
2279 csa_history=prefix(:lenpre)//'.CSA.history'
2280 csa_bank=prefix(:lenpre)//'.CSA.bank'
2281 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2282 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2283 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2284 c!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2285 csa_int=prefix(:lenpre)//'.int'
2286 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2287 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2288 csa_in=prefix(:lenpre)//'.CSA.in'
2289 c print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2292 write (iout,'(80(1h-))')
2293 write (iout,'(30x,a)') "FILE ASSIGNMENT"
2294 write (iout,'(80(1h-))')
2295 write (iout,*) "Input file : ",
2296 & pref_orig(:ilen(pref_orig))//'.inp'
2297 write (iout,*) "Output file : ",
2298 & outname(:ilen(outname))
2300 write (iout,*) "Sidechain potential file : ",
2301 & sidename(:ilen(sidename))
2303 write (iout,*) "SCp potential file : ",
2304 & scpname(:ilen(scpname))
2306 write (iout,*) "Electrostatic potential file : ",
2307 & elename(:ilen(elename))
2308 write (iout,*) "Cumulant coefficient file : ",
2309 & fouriername(:ilen(fouriername))
2310 write (iout,*) "Torsional parameter file : ",
2311 & torname(:ilen(torname))
2312 write (iout,*) "Double torsional parameter file : ",
2313 & tordname(:ilen(tordname))
2314 write (iout,*) "SCCOR parameter file : ",
2315 & sccorname(:ilen(sccorname))
2316 write (iout,*) "Bond & inertia constant file : ",
2317 & bondname(:ilen(bondname))
2318 write (iout,*) "Bending parameter file : ",
2319 & thetname(:ilen(thetname))
2320 write (iout,*) "Rotamer parameter file : ",
2321 & rotname(:ilen(rotname))
2322 write (iout,*) "Thetpdb parameter file : ",
2323 & thetname_pdb(:ilen(thetname_pdb))
2324 write (iout,*) "Threading database : ",
2325 & patname(:ilen(patname))
2327 &write (iout,*)" DIRTMP : ",
2329 write (iout,'(80(1h-))')
2333 c----------------------------------------------------------------------------
2334 subroutine card_concat(card)
2335 implicit real*8 (a-h,o-z)
2336 include 'DIMENSIONS'
2337 include 'COMMON.IOUNITS'
2339 character*80 karta,ucase
2341 read (inp,'(a)') karta
2344 do while (karta(80:80).eq.'&')
2345 card=card(:ilen(card)+1)//karta(:79)
2346 read (inp,'(a)') karta
2349 card=card(:ilen(card)+1)//karta
2352 c----------------------------------------------------------------------------------
2354 implicit real*8 (a-h,o-z)
2355 include 'DIMENSIONS'
2356 include 'COMMON.CHAIN'
2357 include 'COMMON.IOUNITS'
2359 include 'COMMON.CONTROL'
2360 open(irest2,file=rest2name,status='unknown')
2361 read(irest2,*) totT,EK,potE,totE,t_bath
2364 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2367 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2369 if(usampl.or.homol_nset.gt.1) then
2370 read (irest2,*) iset
2375 c---------------------------------------------------------------------------------
2376 subroutine read_fragments
2377 implicit real*8 (a-h,o-z)
2378 include 'DIMENSIONS'
2382 include 'COMMON.SETUP'
2383 include 'COMMON.CHAIN'
2384 include 'COMMON.IOUNITS'
2386 include 'COMMON.CONTROL'
2388 read(inp,*) nset,nfrag,npair,nfrag_back
2389 if(me.eq.king.or..not.out1file)
2390 & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2391 & " nfrag_back",nfrag_back
2393 read(inp,*) mset(iset1)
2395 read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2397 if(me.eq.king.or..not.out1file)
2398 & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2399 & ifrag(2,i,iset1), qinfrag(i,iset1)
2402 read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2404 if(me.eq.king.or..not.out1file)
2405 & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2406 & ipair(2,i,iset1), qinpair(i,iset1)
2409 read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2410 & wfrag_back(3,i,iset1),
2411 & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2412 if(me.eq.king.or..not.out1file)
2413 & write(iout,*) "A",i,wfrag_back(1,i,iset1),
2414 & wfrag_back(2,i,iset1),
2415 & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2416 & ifrag_back(2,i,iset1)
2421 C---------------------------------------------------------------------------
2422 subroutine read_afminp
2423 implicit real*8 (a-h,o-z)
2424 include 'DIMENSIONS'
2428 include 'COMMON.SETUP'
2429 include 'COMMON.CONTROL'
2430 include 'COMMON.CHAIN'
2431 include 'COMMON.IOUNITS'
2432 include 'COMMON.SBRIDGE'
2433 character*320 afmcard
2434 c print *, "wchodze"
2435 call card_concat(afmcard)
2436 call readi(afmcard,"BEG",afmbeg,0)
2437 call readi(afmcard,"END",afmend,0)
2438 call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
2439 call reada(afmcard,"VEL",velAFMconst,0.0d0)
2440 print *,'FORCE=' ,forceAFMconst
2441 CCCC NOW PROPERTIES FOR AFM
2444 distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
2446 distafminit=dsqrt(distafminit)
2447 print *,'initdist',distafminit
2451 c-------------------------------------------------------------------------------
2452 subroutine read_dist_constr
2453 implicit real*8 (a-h,o-z)
2454 include 'DIMENSIONS'
2458 include 'COMMON.SETUP'
2459 include 'COMMON.CONTROL'
2460 include 'COMMON.CHAIN'
2461 include 'COMMON.IOUNITS'
2462 include 'COMMON.SBRIDGE'
2463 integer ifrag_(2,100),ipair_(2,100)
2464 double precision wfrag_(100),wpair_(100)
2465 character*500 controlcard
2466 c write (iout,*) "Calling read_dist_constr"
2467 c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2469 call card_concat(controlcard)
2470 call readi(controlcard,"NFRAG",nfrag_,0)
2471 call readi(controlcard,"NPAIR",npair_,0)
2472 call readi(controlcard,"NDIST",ndist_,0)
2473 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2474 call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2475 call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2476 call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2477 call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2478 c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2479 c write (iout,*) "IFRAG"
2481 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2483 c write (iout,*) "IPAIR"
2485 c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2487 if (.not.refstr .and. nfrag.gt.0) then
2489 & "ERROR: no reference structure to compute distance restraints"
2491 & "Restraints must be specified explicitly (NDIST=number)"
2494 if (nfrag.lt.2 .and. npair.gt.0) then
2495 write (iout,*) "ERROR: Less than 2 fragments specified",
2496 & " but distance restraints between pairs requested"
2501 if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2502 if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2503 & ifrag_(2,i)=nstart_sup+nsup-1
2504 c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2506 if (wfrag_(i).gt.0.0d0) then
2507 do j=ifrag_(1,i),ifrag_(2,i)-1
2508 do k=j+1,ifrag_(2,i)
2509 c write (iout,*) "j",j," k",k
2511 if (constr_dist.eq.1) then
2516 forcon(nhpb)=wfrag_(i)
2517 else if (constr_dist.eq.2) then
2518 if (ddjk.le.dist_cut) then
2523 forcon(nhpb)=wfrag_(i)
2530 forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2533 if (.not.out1file .or. me.eq.king)
2534 & write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2535 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2537 write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2538 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2545 if (wpair_(i).gt.0.0d0) then
2553 do j=ifrag_(1,ii),ifrag_(2,ii)
2554 do k=ifrag_(1,jj),ifrag_(2,jj)
2558 forcon(nhpb)=wpair_(i)
2559 dhpb(nhpb)=dist(j,k)
2561 if (.not.out1file .or. me.eq.king)
2562 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2563 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2565 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2566 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2573 read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2574 if (forcon(nhpb+1).gt.0.0d0) then
2576 dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2578 if (.not.out1file .or. me.eq.king)
2579 & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2580 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2582 write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2583 & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2590 c-------------------------------------------------------------------------------
2592 subroutine read_constr_homology
2594 include 'DIMENSIONS'
2598 include 'COMMON.SETUP'
2599 include 'COMMON.CONTROL'
2600 include 'COMMON.CHAIN'
2601 include 'COMMON.IOUNITS'
2603 include 'COMMON.GEO'
2604 include 'COMMON.INTERACT'
2605 include 'COMMON.NAMES'
2607 c For new homol impl
2609 include 'COMMON.VAR'
2612 c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2614 c common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2615 c & sigma_odl_temp(maxres,maxres,max_template)
2617 character*24 model_ki_dist, model_ki_angle
2618 character*500 controlcard
2619 integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2620 logical lprn /.true./
2625 c FP - Nov. 2014 Temporary specifications for new vars
2627 double precision rescore_tmp,x12,y12,z12,rescore2_tmp,
2629 double precision, dimension (max_template,maxres) :: rescore
2630 double precision, dimension (max_template,maxres) :: rescore2
2631 double precision, dimension (max_template,maxres) :: rescore3
2632 character*24 pdbfile,tpl_k_rescore
2633 c -----------------------------------------------------------------
2634 c Reading multiple PDB ref structures and calculation of retraints
2635 c not using pre-computed ones stored in files model_ki_{dist,angle}
2637 c -----------------------------------------------------------------
2640 c Alternative: reading from input
2641 call card_concat(controlcard)
2642 call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2643 call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2644 call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2645 call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2646 call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2647 call reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0)
2648 call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2649 read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2650 start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0)
2651 if(.not.read2sigma.and.start_from_model) then
2652 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
2653 & write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA'
2654 start_from_model=.false.
2656 if(start_from_model .and. (me.eq.king .or. .not. out1file))
2657 & write(iout,*) 'START_FROM_MODELS is ON'
2658 if(start_from_model .and. rest) then
2659 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2660 write(iout,*) 'START_FROM_MODELS is OFF'
2661 write(iout,*) 'remove restart keyword from input'
2664 if (homol_nset.gt.1)then
2665 call card_concat(controlcard)
2666 read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2667 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2668 write(iout,*) "iset homology_weight "
2670 write(iout,*) i,waga_homology(i)
2673 iset=mod(kolor,homol_nset)+1
2676 waga_homology(1)=1.0
2679 cd write (iout,*) "nnt",nnt," nct",nct
2686 c write(iout,*) 'nnt=',nnt,'nct=',nct
2689 do k=1,constr_homology
2702 do k=1,constr_homology
2704 read(inp,'(a)') pdbfile
2705 if(me.eq.king .or. .not. out1file)
2706 & write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2707 & pdbfile(:ilen(pdbfile))
2708 open(ipdbin,file=pdbfile,status='old',err=33)
2710 33 write (iout,'(a,5x,a)') 'Error opening PDB file',
2711 & pdbfile(:ilen(pdbfile))
2714 c print *,'Begin reading pdb data'
2716 c Files containing res sim or local scores (former containing sigmas)
2719 write(kic2,'(bz,i2.2)') k
2721 tpl_k_rescore="template"//kic2//".sco"
2724 if (read2sigma) then
2725 call readpdb_template(k)
2730 c Distance restraints
2733 C Copy the coordinates from reference coordinates (?)
2737 c write (iout,*) "c(",j,i,") =",c(j,i)
2741 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2743 c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2744 open (ientin,file=tpl_k_rescore,status='old')
2745 if (nnt.gt.1) rescore(k,1)=0.0d0
2746 do irec=nnt,nct ! loop for reading res sim
2747 if (read2sigma) then
2748 read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2749 & rescore3_tmp,idomain_tmp
2751 idomain(k,i_tmp)=idomain_tmp
2752 rescore(k,i_tmp)=rescore_tmp
2753 rescore2(k,i_tmp)=rescore2_tmp
2754 rescore3(k,i_tmp)=rescore3_tmp
2755 if (.not. out1file .or. me.eq.king)
2756 & write(iout,'(a7,i5,3f10.5,i5)') "rescore",
2757 & i_tmp,rescore2_tmp,rescore_tmp,
2758 & rescore3_tmp,idomain_tmp
2761 read (ientin,*,end=1401) rescore_tmp
2763 c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2764 rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2765 c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2770 if (waga_dist.ne.0.0d0) then
2778 distal=dsqrt(x12*x12+y12*y12+z12*z12)
2779 c write (iout,*) k,i,j,distal,dist2_cut
2781 if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0
2782 & .and. distal.le.dist2_cut ) then
2788 c write (iout,*) "k",k
2789 c write (iout,*) "i",i," j",j," constr_homology",
2794 if (read2sigma) then
2797 sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2799 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2800 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2801 & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2803 if (odl(k,ii).le.dist_cut) then
2804 sigma_odl(k,ii)=rescore(k,i)+rescore(k,j)
2807 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2808 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2810 sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2811 & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2815 sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii))
2818 l_homo(k,ii)=.false.
2825 c Theta, dihedral and SC retraints
2827 if (waga_angle.gt.0.0d0) then
2828 c open (ientin,file=tpl_k_sigma_dih,status='old')
2829 c do irec=1,maxres-3 ! loop for reading sigma_dih
2830 c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2831 c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2832 c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2833 c & sigma_dih(k,i+nnt-1)
2838 if (idomain(k,i).eq.0) then
2842 dih(k,i)=phiref(i) ! right?
2843 c read (ientin,*) sigma_dih(k,i) ! original variant
2844 c write (iout,*) "dih(",k,i,") =",dih(k,i)
2845 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2846 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2847 c & "rescore(",k,i-2,") =",rescore(k,i-2),
2848 c & "rescore(",k,i-3,") =",rescore(k,i-3)
2850 sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2851 & rescore(k,i-2)+rescore(k,i-3))/4.0
2852 c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2853 c write (iout,*) "Raw sigmas for dihedral angle restraints"
2854 c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2855 c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2856 c rescore(k,i-2)*rescore(k,i-3) ! right expression ?
2857 c Instead of res sim other local measure of b/b str reliability possible
2858 if (sigma_dih(k,i).ne.0)
2859 & sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2860 c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2865 if (waga_theta.gt.0.0d0) then
2866 c open (ientin,file=tpl_k_sigma_theta,status='old')
2867 c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2868 c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2869 c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2870 c & sigma_theta(k,i+nnt-1)
2875 do i = nnt+2,nct ! right? without parallel.
2876 c do i = i=1,nres ! alternative for bounds acc to readpdb?
2877 c do i=ithet_start,ithet_end ! with FG parallel.
2878 if (idomain(k,i).eq.0) then
2879 sigma_theta(k,i)=0.0
2882 thetatpl(k,i)=thetaref(i)
2883 c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2884 c write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2885 c & "rescore(",k,i-1,") =",rescore(k,i-1),
2886 c & "rescore(",k,i-2,") =",rescore(k,i-2)
2887 c read (ientin,*) sigma_theta(k,i) ! 1st variant
2888 sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2889 & rescore(k,i-2))/3.0
2890 c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2891 if (sigma_theta(k,i).ne.0)
2892 & sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2894 c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2895 c rescore(k,i-2) ! right expression ?
2896 c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2900 if (waga_d.gt.0.0d0) then
2901 c open (ientin,file=tpl_k_sigma_d,status='old')
2902 c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2903 c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2904 c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2905 c & sigma_d(k,i+nnt-1)
2909 do i = nnt,nct ! right? without parallel.
2910 c do i=2,nres-1 ! alternative for bounds acc to readpdb?
2911 c do i=loc_start,loc_end ! with FG parallel.
2912 if (itype(i).eq.10) cycle
2913 if (idomain(k,i).eq.0 ) then
2920 c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2921 c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2922 c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2923 c write(iout,*) "rescore(",k,i,") =",rescore(k,i)
2924 sigma_d(k,i)=rescore3(k,i) ! right expression ?
2925 if (sigma_d(k,i).ne.0)
2926 & sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2928 c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ?
2929 c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2930 c read (ientin,*) sigma_d(k,i) ! 1st variant
2935 c remove distance restraints not used in any model from the list
2936 c shift data in all arrays
2938 if (waga_dist.ne.0.0d0) then
2944 if (ii_in_use(ii).eq.0.and.liiflag) then
2948 if (ii_in_use(ii).ne.0.and..not.liiflag.or.
2949 & .not.liiflag.and.ii.eq.lim_odl) then
2950 if (ii.eq.lim_odl) then
2951 iishift=ii-iistart+1
2956 do ki=iistart,lim_odl-iishift
2957 ires_homo(ki)=ires_homo(ki+iishift)
2958 jres_homo(ki)=jres_homo(ki+iishift)
2959 ii_in_use(ki)=ii_in_use(ki+iishift)
2960 do k=1,constr_homology
2961 odl(k,ki)=odl(k,ki+iishift)
2962 sigma_odl(k,ki)=sigma_odl(k,ki+iishift)
2963 l_homo(k,ki)=l_homo(k,ki+iishift)
2967 lim_odl=lim_odl-iishift
2972 if (constr_homology.gt.0) call homology_partition
2973 if (constr_homology.gt.0) call init_int_table
2974 c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2975 c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2979 if (.not.lprn) return
2980 cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2981 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2982 write (iout,*) "Distance restraints from templates"
2984 write(iout,'(3i5,100(2f8.2,1x,l1,4x))')
2985 & ii,ires_homo(ii),jres_homo(ii),
2986 & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
2987 & ki=1,constr_homology)
2989 write (iout,*) "Dihedral angle restraints from templates"
2991 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
2992 & (rad2deg*dih(ki,i),
2993 & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2995 write (iout,*) "Virtual-bond angle restraints from templates"
2997 write (iout,'(i5,a4,100(2f8.2,4x))') i,restyp(itype(i)),
2998 & (rad2deg*thetatpl(ki,i),
2999 & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3001 write (iout,*) "SC restraints from templates"
3003 write(iout,'(i5,100(4f8.2,4x))') i,
3004 & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3005 & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3008 c -----------------------------------------------------------------
3011 c----------------------------------------------------------------------
3014 subroutine flush(iu)
3019 subroutine flush(iu)
3024 c------------------------------------------------------------------------------
3025 subroutine copy_to_tmp(source)
3026 include "DIMENSIONS"
3027 include "COMMON.IOUNITS"
3028 character*(*) source
3029 character* 256 tmpfile
3033 tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3034 inquire(file=tmpfile,exist=ex)
3036 write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3037 & " to temporary directory..."
3038 write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3039 call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3043 c------------------------------------------------------------------------------
3044 subroutine move_from_tmp(source)
3045 include "DIMENSIONS"
3046 include "COMMON.IOUNITS"
3047 character*(*) source
3050 write (*,*) "Moving ",source(:ilen(source)),
3051 & " from temporary directory to working directory"
3052 write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3053 call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3056 c------------------------------------------------------------------------------
3057 subroutine random_init(seed)
3059 C Initialize random number generator
3061 implicit real*8 (a-h,o-z)
3062 include 'DIMENSIONS'
3068 logical OKRandom, prng_restart
3070 integer iseed_array(4)
3072 include 'COMMON.IOUNITS'
3073 include 'COMMON.TIME1'
3074 include 'COMMON.THREAD'
3075 include 'COMMON.SBRIDGE'
3076 include 'COMMON.CONTROL'
3077 include 'COMMON.MCM'
3078 include 'COMMON.MAP'
3079 include 'COMMON.HEADER'
3080 include 'COMMON.CSA'
3081 include 'COMMON.CHAIN'
3082 include 'COMMON.MUCA'
3084 include 'COMMON.FFIELD'
3085 include 'COMMON.SETUP'
3086 iseed=-dint(dabs(seed))
3087 if (iseed.eq.0) then
3088 write (iout,'(/80(1h*)/20x,a/80(1h*))')
3089 & 'Random seed undefined. The program will stop.'
3090 write (*,'(/80(1h*)/20x,a/80(1h*))')
3091 & 'Random seed undefined. The program will stop.'
3093 call mpi_finalize(mpi_comm_world,ierr)
3095 stop 'Bad random seed.'
3098 if (fg_rank.eq.0) then
3102 if(me.eq.king .or. .not. out1file)
3103 & write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3104 write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3105 OKRandom = prng_restart(me,iseedi8)
3108 tmp=65536.0d0**(4-i)
3109 iseed_array(i) = dint(seed/tmp)
3110 seed=seed-iseed_array(i)*tmp
3113 & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3114 & (iseed_array(i),i=1,4)
3115 write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3116 & (iseed_array(i),i=1,4)
3117 OKRandom = prng_restart(me,iseed_array)
3120 c r1 = prng_next(me)
3121 r1=ran_number(0.0D0,1.0D0)
3122 if(me.eq.king .or. .not. out1file)
3123 & write (iout,*) 'ran_num',r1
3124 if (r1.lt.0.0d0) OKRandom=.false.
3126 if (.not.OKRandom) then
3127 write (iout,*) 'PRNG IS NOT WORKING!!!'
3128 print *,'PRNG IS NOT WORKING!!!'
3131 call mpi_abort(mpi_comm_world,error_msg,ierr)
3134 write (iout,*) 'too many processors for parallel prng'
3135 write (*,*) 'too many processors for parallel prng'
3143 write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)