8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm
743 integer :: maxinter,junk,kk,ii,ncatprotparm
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2
750 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
751 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
752 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 ! For printing parameters after they are read set the following in the UNRES
757 ! setenv PRINT_PARM YES
759 ! To print parameters in LaTeX format rather than as ASCII tables:
763 call getenv_loc("PRINT_PARM",lancuch)
764 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
765 call getenv_loc("LATEX",lancuch)
766 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
768 dwa16=2.0d0**(1.0d0/6.0d0)
770 ! Assign virtual-bond length
773 vblinv2=vblinv*vblinv
775 ! Read the virtual-bond parameters, masses, and moments of inertia
776 ! and Stokes' radii of the peptide group and side chains
778 allocate(dsc(ntyp1)) !(ntyp1)
779 allocate(dsc_inv(ntyp1)) !(ntyp1)
780 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
781 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
782 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
783 allocate(nbondterm(ntyp)) !(ntyp)
784 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
785 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
786 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(long_r_sidechain(ntyp))
788 allocate(short_r_sidechain(ntyp))
793 allocate(msc(ntyp+1)) !(ntyp+1)
794 allocate(isc(ntyp+1)) !(ntyp+1)
795 allocate(restok(ntyp+1)) !(ntyp+1)
797 read (ibond,*) vbldp0,akp,mp,ip,pstok
800 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
801 dsc(i) = vbldsc0(1,i)
805 dsc_inv(i)=1.0D0/dsc(i)
814 allocate(msc(ntyp+1,5)) !(ntyp+1)
815 allocate(isc(ntyp+1,5)) !(ntyp+1)
816 allocate(restok(ntyp+1,5)) !(ntyp+1)
818 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
820 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
821 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
822 dsc(i) = vbldsc0(1,i)
826 dsc_inv(i)=1.0D0/dsc(i)
830 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
833 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
834 ! dsc(i) = vbldsc0_nucl(1,i)
838 ! dsc_inv(i)=1.0D0/dsc(i)
841 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
842 ! do i=1,ntyp_molec(2)
843 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
844 ! aksc_nucl(j,i),abond0_nucl(j,i),&
845 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
846 ! dsc(i) = vbldsc0(1,i)
850 ! dsc_inv(i)=1.0D0/dsc(i)
855 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
856 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
858 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
860 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
861 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
863 write (iout,'(13x,3f10.5)') &
864 vbldsc0(j,i),aksc(j,i),abond0(j,i)
869 read(iion,*) msc(i,5),restok(i,5)
870 print *,msc(i,5),restok(i,5)
874 read (iion,*) ncatprotparm
875 allocate(catprm(ncatprotparm,4))
877 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
880 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
881 !----------------------------------------------------
882 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
883 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
884 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
885 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
886 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
887 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
897 allocate(liptranene(ntyp))
898 !C reading lipid parameters
899 write (iout,*) "iliptranpar",iliptranpar
901 read(iliptranpar,*) pepliptran
904 read(iliptranpar,*) liptranene(i)
905 print *,liptranene(i)
911 ! Read the parameters of the probability distribution/energy expression
912 ! of the virtual-bond valence angles theta
915 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
916 (bthet(j,i,1,1),j=1,2)
917 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
918 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
919 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
923 athet(1,i,1,-1)=athet(1,i,1,1)
924 athet(2,i,1,-1)=athet(2,i,1,1)
925 bthet(1,i,1,-1)=-bthet(1,i,1,1)
926 bthet(2,i,1,-1)=-bthet(2,i,1,1)
927 athet(1,i,-1,1)=-athet(1,i,1,1)
928 athet(2,i,-1,1)=-athet(2,i,1,1)
929 bthet(1,i,-1,1)=bthet(1,i,1,1)
930 bthet(2,i,-1,1)=bthet(2,i,1,1)
934 athet(1,i,-1,-1)=athet(1,-i,1,1)
935 athet(2,i,-1,-1)=-athet(2,-i,1,1)
936 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
937 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
938 athet(1,i,-1,1)=athet(1,-i,1,1)
939 athet(2,i,-1,1)=-athet(2,-i,1,1)
940 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
941 bthet(2,i,-1,1)=bthet(2,-i,1,1)
942 athet(1,i,1,-1)=-athet(1,-i,1,1)
943 athet(2,i,1,-1)=athet(2,-i,1,1)
944 bthet(1,i,1,-1)=bthet(1,-i,1,1)
945 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
950 polthet(j,i)=polthet(j,-i)
953 gthet(j,i)=gthet(j,-i)
961 'Parameters of the virtual-bond valence angles:'
962 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
963 ' ATHETA0 ',' A1 ',' A2 ',&
966 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
967 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
969 write (iout,'(/a/9x,5a/79(1h-))') &
970 'Parameters of the expression for sigma(theta_c):',&
971 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
972 ' ALPH3 ',' SIGMA0C '
974 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
975 (polthet(j,i),j=0,3),sigc0(i)
977 write (iout,'(/a/9x,5a/79(1h-))') &
978 'Parameters of the second gaussian:',&
979 ' THETA0 ',' SIGMA0 ',' G1 ',&
982 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
983 sig0(i),(gthet(j,i),j=1,3)
987 'Parameters of the virtual-bond valence angles:'
988 write (iout,'(/a/9x,5a/79(1h-))') &
989 'Coefficients of expansion',&
990 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
991 ' b1*10^1 ',' b2*10^1 '
993 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
994 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
995 (10*bthet(j,i,1,1),j=1,2)
997 write (iout,'(/a/9x,5a/79(1h-))') &
998 'Parameters of the expression for sigma(theta_c):',&
999 ' alpha0 ',' alph1 ',' alph2 ',&
1000 ' alhp3 ',' sigma0c '
1002 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1003 (polthet(j,i),j=0,3),sigc0(i)
1005 write (iout,'(/a/9x,5a/79(1h-))') &
1006 'Parameters of the second gaussian:',&
1007 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1010 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1011 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1017 ! Read the parameters of Utheta determined from ab initio surfaces
1018 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1020 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1021 ntheterm3,nsingle,ndouble
1022 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1024 !----------------------------------------------------
1025 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1026 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1027 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1028 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1029 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1032 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1033 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1034 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1035 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1036 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1037 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1038 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1039 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1040 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1041 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1042 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1043 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1044 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1045 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1046 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1047 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1048 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1050 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1052 ithetyp(i)=-ithetyp(-i)
1055 aa0thet(:,:,:,:)=0.0d0
1056 aathet(:,:,:,:,:)=0.0d0
1057 bbthet(:,:,:,:,:,:)=0.0d0
1058 ccthet(:,:,:,:,:,:)=0.0d0
1059 ddthet(:,:,:,:,:,:)=0.0d0
1060 eethet(:,:,:,:,:,:)=0.0d0
1061 ffthet(:,:,:,:,:,:,:)=0.0d0
1062 ggthet(:,:,:,:,:,:,:)=0.0d0
1064 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1066 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1067 ! VAR:1=non-glicyne non-proline 2=proline
1068 ! VAR:negative values for D-aminoacid
1070 do j=-nthetyp,nthetyp
1071 do k=-nthetyp,nthetyp
1072 read (ithep,'(6a)',end=111,err=111) res1
1073 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1074 ! VAR: aa0thet is variable describing the average value of Foureir
1075 ! VAR: expansion series
1076 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1077 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1078 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1079 read (ithep,*,end=111,err=111) &
1080 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1081 read (ithep,*,end=111,err=111) &
1082 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1083 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1084 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1085 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1087 read (ithep,*,end=111,err=111) &
1088 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1089 ffthet(lll,llll,ll,i,j,k,iblock),&
1090 ggthet(llll,lll,ll,i,j,k,iblock),&
1091 ggthet(lll,llll,ll,i,j,k,iblock),&
1092 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1097 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1098 ! coefficients of theta-and-gamma-dependent terms are zero.
1099 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1100 ! RECOMENTDED AFTER VERSION 3.3)
1104 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1105 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1107 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1108 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1111 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1113 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1116 ! AND COMMENT THE LOOPS BELOW
1120 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1121 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1123 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1124 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1127 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1129 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1134 ! Substitution for D aminoacids from symmetry.
1137 do j=-nthetyp,nthetyp
1138 do k=-nthetyp,nthetyp
1139 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1141 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1145 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1146 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1147 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1148 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1154 ffthet(llll,lll,ll,i,j,k,iblock)= &
1155 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1156 ffthet(lll,llll,ll,i,j,k,iblock)= &
1157 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1158 ggthet(llll,lll,ll,i,j,k,iblock)= &
1159 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1160 ggthet(lll,llll,ll,i,j,k,iblock)= &
1161 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1170 ! Control printout of the coefficients of virtual-bond-angle potentials
1173 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1178 write (iout,'(//4a)') &
1179 'Type ',onelett(i),onelett(j),onelett(k)
1180 write (iout,'(//a,10x,a)') " l","a[l]"
1181 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1182 write (iout,'(i2,1pe15.5)') &
1183 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1185 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1186 "b",l,"c",l,"d",l,"e",l
1188 write (iout,'(i2,4(1pe15.5))') m,&
1189 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1190 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1194 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1195 "f+",l,"f-",l,"g+",l,"g-",l
1198 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1199 ffthet(n,m,l,i,j,k,iblock),&
1200 ffthet(m,n,l,i,j,k,iblock),&
1201 ggthet(n,m,l,i,j,k,iblock),&
1202 ggthet(m,n,l,i,j,k,iblock)
1212 write (2,*) "Start reading THETA_PDB",ithep_pdb
1214 ! write (2,*) 'i=',i
1215 read (ithep_pdb,*,err=111,end=111) &
1216 a0thet(i),(athet(j,i,1,1),j=1,2),&
1217 (bthet(j,i,1,1),j=1,2)
1218 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1219 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1220 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1221 sigc0(i)=sigc0(i)**2
1224 athet(1,i,1,-1)=athet(1,i,1,1)
1225 athet(2,i,1,-1)=athet(2,i,1,1)
1226 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1227 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1228 athet(1,i,-1,1)=-athet(1,i,1,1)
1229 athet(2,i,-1,1)=-athet(2,i,1,1)
1230 bthet(1,i,-1,1)=bthet(1,i,1,1)
1231 bthet(2,i,-1,1)=bthet(2,i,1,1)
1234 a0thet(i)=a0thet(-i)
1235 athet(1,i,-1,-1)=athet(1,-i,1,1)
1236 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1237 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1238 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1239 athet(1,i,-1,1)=athet(1,-i,1,1)
1240 athet(2,i,-1,1)=-athet(2,-i,1,1)
1241 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1242 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1243 athet(1,i,1,-1)=-athet(1,-i,1,1)
1244 athet(2,i,1,-1)=athet(2,-i,1,1)
1245 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1246 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1247 theta0(i)=theta0(-i)
1251 polthet(j,i)=polthet(j,-i)
1254 gthet(j,i)=gthet(j,-i)
1257 write (2,*) "End reading THETA_PDB"
1261 !--------------- Reading theta parameters for nucleic acid-------
1262 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1263 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1264 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1265 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1266 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1267 nthetyp_nucl+1,nthetyp_nucl+1))
1268 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1269 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1270 nthetyp_nucl+1,nthetyp_nucl+1))
1271 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1272 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1273 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1274 nthetyp_nucl+1,nthetyp_nucl+1))
1275 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1276 nthetyp_nucl+1,nthetyp_nucl+1))
1277 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1278 nthetyp_nucl+1,nthetyp_nucl+1))
1279 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1280 nthetyp_nucl+1,nthetyp_nucl+1))
1281 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1282 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1283 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1284 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1285 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1286 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1288 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1289 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1291 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1293 aa0thet_nucl(:,:,:)=0.0d0
1294 aathet_nucl(:,:,:,:)=0.0d0
1295 bbthet_nucl(:,:,:,:,:)=0.0d0
1296 ccthet_nucl(:,:,:,:,:)=0.0d0
1297 ddthet_nucl(:,:,:,:,:)=0.0d0
1298 eethet_nucl(:,:,:,:,:)=0.0d0
1299 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1300 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1305 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1306 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1307 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1308 read (ithep_nucl,*,end=111,err=111) &
1309 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1310 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1311 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1312 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1313 read (ithep_nucl,*,end=111,err=111) &
1314 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1315 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1316 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1321 !-------------------------------------------
1322 allocate(nlob(ntyp1)) !(ntyp1)
1323 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1324 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1325 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1332 gaussc(:,:,:,:)=0.0D0
1336 ! Read the parameters of the probability distribution/energy expression
1337 ! of the side chains.
1340 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1344 dsc_inv(i)=1.0D0/dsc(i)
1355 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1356 ((blower(k,l,1),l=1,k),k=1,3)
1357 censc(1,1,-i)=censc(1,1,i)
1358 censc(2,1,-i)=censc(2,1,i)
1359 censc(3,1,-i)=-censc(3,1,i)
1361 read (irotam,*,end=112,err=112) bsc(j,i)
1362 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1363 ((blower(k,l,j),l=1,k),k=1,3)
1364 censc(1,j,-i)=censc(1,j,i)
1365 censc(2,j,-i)=censc(2,j,i)
1366 censc(3,j,-i)=-censc(3,j,i)
1367 ! BSC is amplitude of Gaussian
1374 akl=akl+blower(k,m,j)*blower(l,m,j)
1378 if (((k.eq.3).and.(l.ne.3)) &
1379 .or.((l.eq.3).and.(k.ne.3))) then
1380 gaussc(k,l,j,-i)=-akl
1381 gaussc(l,k,j,-i)=-akl
1383 gaussc(k,l,j,-i)=akl
1384 gaussc(l,k,j,-i)=akl
1393 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1396 if (nlobi.gt.0) then
1398 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1399 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1400 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1401 'log h',(bsc(j,i),j=1,nlobi)
1402 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1403 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1405 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1406 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1409 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1410 write (iout,'(a,f10.4,4(16x,f10.4))') &
1411 'Center ',(bsc(j,i),j=1,nlobi)
1412 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1421 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1422 ! added by Urszula Kozlowska 07/11/2007
1424 !el Maximum number of SC local term fitting function coefficiants
1425 !el integer,parameter :: maxsccoef=65
1427 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1430 read (irotam,*,end=112,err=112)
1432 read (irotam,*,end=112,err=112)
1435 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1439 !---------reading nucleic acid parameters for rotamers-------------------
1440 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1441 do i=1,ntyp_molec(2)
1442 read (irotam_nucl,*,end=112,err=112)
1444 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1450 write (iout,*) "Base rotamer parameters"
1451 do i=1,ntyp_molec(2)
1452 write (iout,'(a)') restyp(i,2)
1453 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1458 ! Read the parameters of the probability distribution/energy expression
1459 ! of the side chains.
1461 write (2,*) "Start reading ROTAM_PDB"
1463 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1467 dsc_inv(i)=1.0D0/dsc(i)
1478 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1479 ((blower(k,l,1),l=1,k),k=1,3)
1481 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1482 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1483 ((blower(k,l,j),l=1,k),k=1,3)
1490 akl=akl+blower(k,m,j)*blower(l,m,j)
1500 write (2,*) "End reading ROTAM_PDB"
1506 ! Read torsional parameters in old format
1508 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1510 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1511 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1512 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1514 !el from energy module--------
1515 allocate(v1(nterm_old,ntortyp,ntortyp))
1516 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1517 !el---------------------------
1522 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1528 write (iout,'(/a/)') 'Torsional constants:'
1531 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1532 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
1538 ! Read torsional parameters
1540 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1541 read (itorp,*,end=113,err=113) ntortyp
1542 !el from energy module---------
1543 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1544 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1546 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1547 allocate(vlor2(maxlor,ntortyp,ntortyp))
1548 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
1549 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1551 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
1552 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1553 !el---------------------------
1556 !el---------------------------
1558 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1560 itortyp(i)=-itortyp(-i)
1562 itortyp(ntyp1)=ntortyp
1563 itortyp(-ntyp1)=-ntortyp
1565 write (iout,*) 'ntortyp',ntortyp
1567 do j=-ntortyp+1,ntortyp-1
1568 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
1570 nterm(-i,-j,iblock)=nterm(i,j,iblock)
1571 nlor(-i,-j,iblock)=nlor(i,j,iblock)
1574 do k=1,nterm(i,j,iblock)
1575 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
1577 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
1578 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
1579 v0ij=v0ij+si*v1(k,i,j,iblock)
1581 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
1582 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
1583 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
1585 do k=1,nlor(i,j,iblock)
1586 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
1587 vlor2(k,i,j),vlor3(k,i,j)
1588 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
1591 v0(-i,-j,iblock)=v0ij
1597 write (iout,'(/a/)') 'Torsional constants:'
1599 do i=-ntortyp,ntortyp
1600 do j=-ntortyp,ntortyp
1601 write (iout,*) 'ityp',i,' jtyp',j
1602 write (iout,*) 'Fourier constants'
1603 do k=1,nterm(i,j,iblock)
1604 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
1607 write (iout,*) 'Lorenz constants'
1608 do k=1,nlor(i,j,iblock)
1609 write (iout,'(3(1pe15.5))') &
1610 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1616 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
1618 ! 6/23/01 Read parameters for double torsionals
1620 !el from energy module------------
1621 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1622 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1623 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1624 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1625 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1626 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1627 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1628 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
1629 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
1630 !---------------------------------
1634 do j=-ntortyp+1,ntortyp-1
1635 do k=-ntortyp+1,ntortyp-1
1636 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
1637 ! write (iout,*) "OK onelett",
1640 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
1641 .or. t3.ne.toronelet(k)) then
1642 write (iout,*) "Error in double torsional parameter file",&
1645 call MPI_Finalize(Ierror)
1647 stop "Error in double torsional parameter file"
1649 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
1650 ntermd_2(i,j,k,iblock)
1651 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
1652 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
1653 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
1654 ntermd_1(i,j,k,iblock))
1655 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
1656 ntermd_1(i,j,k,iblock))
1657 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
1658 ntermd_1(i,j,k,iblock))
1659 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
1660 ntermd_1(i,j,k,iblock))
1661 ! Martix of D parameters for one dimesional foureir series
1662 do l=1,ntermd_1(i,j,k,iblock)
1663 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
1664 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
1665 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
1666 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
1667 ! write(iout,*) "whcodze" ,
1668 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
1670 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
1671 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
1672 v2s(m,l,i,j,k,iblock),&
1673 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
1674 ! Martix of D parameters for two dimesional fourier series
1675 do l=1,ntermd_2(i,j,k,iblock)
1677 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
1678 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
1679 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
1680 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
1689 write (iout,*) 'Constants for double torsionals'
1692 do j=-ntortyp+1,ntortyp-1
1693 do k=-ntortyp+1,ntortyp-1
1694 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
1695 ' nsingle',ntermd_1(i,j,k,iblock),&
1696 ' ndouble',ntermd_2(i,j,k,iblock)
1698 write (iout,*) 'Single angles:'
1699 do l=1,ntermd_1(i,j,k,iblock)
1700 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
1701 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
1702 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
1703 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
1706 write (iout,*) 'Pairs of angles:'
1707 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1708 do l=1,ntermd_2(i,j,k,iblock)
1709 write (iout,'(i5,20f10.5)') &
1710 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
1713 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
1714 do l=1,ntermd_2(i,j,k,iblock)
1715 write (iout,'(i5,20f10.5)') &
1716 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
1717 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
1726 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1727 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
1728 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
1729 !el from energy module---------
1730 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1731 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1733 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
1734 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
1735 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
1736 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
1738 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
1739 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
1740 !el---------------------------
1743 !el--------------------
1744 read (itorp_nucl,*,end=113,err=113) &
1745 (itortyp_nucl(i),i=1,ntyp_molec(2))
1746 ! print *,itortyp_nucl(:)
1747 !c write (iout,*) 'ntortyp',ntortyp
1750 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
1751 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
1754 do k=1,nterm_nucl(i,j)
1755 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
1756 v0ij=v0ij+si*v1_nucl(k,i,j)
1759 do k=1,nlor_nucl(i,j)
1760 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
1761 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
1762 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
1768 ! Read of Side-chain backbone correlation parameters
1769 ! Modified 11 May 2012 by Adasko
1772 read (isccor,*,end=119,err=119) nsccortyp
1774 !el from module energy-------------
1775 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1776 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
1777 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
1778 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
1779 !-----------------------------------
1781 !el from module energy-------------
1782 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
1784 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1786 isccortyp(i)=-isccortyp(-i)
1788 iscprol=isccortyp(20)
1789 ! write (iout,*) 'ntortyp',ntortyp
1791 !c maxinter is maximum interaction sites
1792 !el from module energy---------
1793 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1794 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1795 -nsccortyp:nsccortyp))
1796 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
1797 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1798 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
1799 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1800 !-----------------------------------
1804 read (isccor,*,end=119,err=119) &
1805 nterm_sccor(i,j),nlor_sccor(i,j)
1811 nterm_sccor(-i,j)=nterm_sccor(i,j)
1812 nterm_sccor(-i,-j)=nterm_sccor(i,j)
1813 nterm_sccor(i,-j)=nterm_sccor(i,j)
1814 do k=1,nterm_sccor(i,j)
1815 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1817 if (j.eq.iscprol) then
1818 if (i.eq.isccortyp(10)) then
1819 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1820 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1822 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
1823 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
1824 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
1825 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
1826 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1827 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1828 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1829 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1832 if (i.eq.isccortyp(10)) then
1833 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
1834 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1836 if (j.eq.isccortyp(10)) then
1837 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
1838 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
1840 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
1841 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
1842 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
1843 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
1844 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
1845 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
1849 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1850 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
1851 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
1852 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
1855 do k=1,nlor_sccor(i,j)
1856 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1857 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1858 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1859 (1+vlor3sccor(k,i,j)**2)
1861 v0sccor(l,i,j)=v0ijsccor
1862 v0sccor(l,-i,j)=v0ijsccor1
1863 v0sccor(l,i,-j)=v0ijsccor2
1864 v0sccor(l,-i,-j)=v0ijsccor3
1870 !el from module energy-------------
1871 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
1873 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
1874 ! write (iout,*) 'ntortyp',ntortyp
1876 !c maxinter is maximum interaction sites
1877 !el from module energy---------
1878 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
1879 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
1880 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
1881 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
1882 !-----------------------------------
1886 read (isccor,*,end=119,err=119) &
1887 nterm_sccor(i,j),nlor_sccor(i,j)
1891 do k=1,nterm_sccor(i,j)
1892 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
1894 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
1897 do k=1,nlor_sccor(i,j)
1898 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
1899 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1900 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
1901 (1+vlor3sccor(k,i,j)**2)
1903 v0sccor(l,i,j)=v0ijsccor !el ,iblock
1912 write (iout,'(/a/)') 'Torsional constants:'
1915 write (iout,*) 'ityp',i,' jtyp',j
1916 write (iout,*) 'Fourier constants'
1917 do k=1,nterm_sccor(i,j)
1918 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
1920 write (iout,*) 'Lorenz constants'
1921 do k=1,nlor_sccor(i,j)
1922 write (iout,'(3(1pe15.5))') &
1923 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
1930 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1931 ! interaction energy of the Gly, Ala, and Pro prototypes.
1936 write (iout,*) "Coefficients of the cumulants"
1938 read (ifourier,*) nloctyp
1939 !write(iout,*) "nloctyp",nloctyp
1940 !el from module energy-------
1941 allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1942 allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
1943 allocate(b1tilde(2,-nloctyp+1:nloctyp+1)) !(2,-maxtor:maxtor)
1944 allocate(cc(2,2,-nloctyp-1:nloctyp+1))
1945 allocate(dd(2,2,-nloctyp-1:nloctyp+1))
1946 allocate(ee(2,2,-nloctyp-1:nloctyp+1))
1947 allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
1948 allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
1952 !--------------------------------
1955 read (ifourier,*,end=115,err=115)
1956 read (ifourier,*,end=115,err=115) (buse(ii),ii=1,13)
1958 write (iout,*) 'Type',i
1959 write (iout,'(a,i2,a,f10.5)') ('buse(',ii,')=',buse(ii),ii=1,13)
1967 B1tilde(1,i) = buse(3)
1968 B1tilde(2,i) =-buse(5)
1969 B1tilde(1,-i) =-buse(3)
1970 B1tilde(2,-i) =buse(5)
1971 ! buse1tilde(1,i)=0.0d0
1972 ! buse1tilde(2,i)=0.0d0
1992 Ctilde(1,1,i)=buse(7)
1993 Ctilde(1,2,i)=buse(9)
1994 Ctilde(2,1,i)=-buse(9)
1995 Ctilde(2,2,i)=buse(7)
1996 Ctilde(1,1,-i)=buse(7)
1997 Ctilde(1,2,-i)=-buse(9)
1998 Ctilde(2,1,-i)=buse(9)
1999 Ctilde(2,2,-i)=buse(7)
2001 ! Ctilde(1,1,i)=0.0d0
2002 ! Ctilde(1,2,i)=0.0d0
2003 ! Ctilde(2,1,i)=0.0d0
2004 ! Ctilde(2,2,i)=0.0d0
2017 Dtilde(1,1,i)=buse(6)
2018 Dtilde(1,2,i)=buse(8)
2019 Dtilde(2,1,i)=-buse(8)
2020 Dtilde(2,2,i)=buse(6)
2021 Dtilde(1,1,-i)=buse(6)
2022 Dtilde(1,2,-i)=-buse(8)
2023 Dtilde(2,1,-i)=buse(8)
2024 Dtilde(2,2,-i)=buse(6)
2026 ! Dtilde(1,1,i)=0.0d0
2027 ! Dtilde(1,2,i)=0.0d0
2028 ! Dtilde(2,1,i)=0.0d0
2029 ! Dtilde(2,2,i)=0.0d0
2030 EE(1,1,i)= buse(10)+buse(11)
2031 EE(2,2,i)=-buse(10)+buse(11)
2032 EE(2,1,i)= buse(12)-buse(13)
2033 EE(1,2,i)= buse(12)+buse(13)
2034 EE(1,1,-i)= buse(10)+buse(11)
2035 EE(2,2,-i)=-buse(10)+buse(11)
2036 EE(2,1,-i)=-buse(12)+buse(13)
2037 EE(1,2,-i)=-buse(12)-buse(13)
2043 ! ee(2,1,i)=ee(1,2,i)
2047 write (iout,*) 'Type',i
2049 write(iout,*) B1(1,i),B1(2,i)
2051 write(iout,*) B2(1,i),B2(2,i)
2054 write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
2058 write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
2062 write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
2067 ! Read electrostatic-interaction parameters
2072 write (iout,'(/a)') 'Electrostatic interaction constants:'
2073 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2074 'IT','JT','APP','BPP','AEL6','AEL3'
2076 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2077 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2078 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2079 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2084 app (i,j)=epp(i,j)*rri*rri
2085 bpp (i,j)=-2.0D0*epp(i,j)*rri
2086 ael6(i,j)=elpp6(i,j)*4.2D0**6
2087 ael3(i,j)=elpp3(i,j)*4.2D0**3
2089 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2095 ! Read side-chain interaction parameters.
2097 !el from module energy - COMMON.INTERACT-------
2098 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2099 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2100 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2102 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2103 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2104 allocate(epslip(ntyp,ntyp))
2112 !--------------------------------
2114 read (isidep,*,end=117,err=117) ipot,expon
2115 if (ipot.lt.1 .or. ipot.gt.5) then
2116 write (iout,'(2a)') 'Error while reading SC interaction',&
2117 'potential file - unknown potential type.'
2119 call MPI_Finalize(Ierror)
2125 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2126 ', exponents are ',expon,2*expon
2127 ! goto (10,20,30,30,40) ipot
2129 !----------------------- LJ potential ---------------------------------
2131 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2132 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2133 (sigma0(i),i=1,ntyp)
2135 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2136 write (iout,'(a/)') 'The epsilon array:'
2137 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2138 write (iout,'(/a)') 'One-body parameters:'
2139 write (iout,'(a,4x,a)') 'residue','sigma'
2140 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2143 !----------------------- LJK potential --------------------------------
2145 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2146 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2147 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2149 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2150 write (iout,'(a/)') 'The epsilon array:'
2151 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2152 write (iout,'(/a)') 'One-body parameters:'
2153 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2154 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2158 !---------------------- GB or BP potential -----------------------------
2161 ! print *,"I AM in SCELE",scelemode
2162 if (scelemode.eq.0) then
2164 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2166 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2167 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2168 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2169 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2171 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2174 ! For the GB potential convert sigma'**2 into chi'
2177 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2181 write (iout,'(/a/)') 'Parameters of the BP potential:'
2182 write (iout,'(a/)') 'The epsilon array:'
2183 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2184 write (iout,'(/a)') 'One-body parameters:'
2185 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2187 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2188 chip(i),alp(i),i=1,ntyp)
2191 ! print *,ntyp,"NTYP"
2192 allocate(icharge(ntyp1))
2193 ! print *,ntyp,icharge(i)
2195 read (isidep,*) (icharge(i),i=1,ntyp)
2196 print *,ntyp,icharge(i)
2197 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2198 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2199 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2200 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2201 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2202 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2203 allocate(epsintab(ntyp,ntyp))
2204 allocate(dtail(2,ntyp,ntyp))
2205 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2206 allocate(wqdip(2,ntyp,ntyp))
2207 allocate(wstate(4,ntyp,ntyp))
2208 allocate(dhead(2,2,ntyp,ntyp))
2209 allocate(nstate(ntyp,ntyp))
2210 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2211 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2214 ! write (*,*) "Im in ALAB", i, " ", j
2216 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2217 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2218 chis(i,j),chis(j,i), &
2219 nstate(i,j),(wstate(k,i,j),k=1,4), &
2220 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2221 dtail(1,i,j),dtail(2,i,j), &
2222 epshead(i,j),sig0head(i,j), &
2223 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2224 alphapol(i,j),alphapol(j,i), &
2225 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j)
2226 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2232 sigma(i,j) = sigma(j,i)
2233 sigmap1(i,j)=sigmap1(j,i)
2234 sigmap2(i,j)=sigmap2(j,i)
2235 sigiso1(i,j)=sigiso1(j,i)
2236 sigiso2(i,j)=sigiso2(j,i)
2237 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2238 nstate(i,j) = nstate(j,i)
2239 dtail(1,i,j) = dtail(1,j,i)
2240 dtail(2,i,j) = dtail(2,j,i)
2242 alphasur(k,i,j) = alphasur(k,j,i)
2243 wstate(k,i,j) = wstate(k,j,i)
2244 alphiso(k,i,j) = alphiso(k,j,i)
2247 dhead(2,1,i,j) = dhead(1,1,j,i)
2248 dhead(2,2,i,j) = dhead(1,2,j,i)
2249 dhead(1,1,i,j) = dhead(2,1,j,i)
2250 dhead(1,2,i,j) = dhead(2,2,j,i)
2252 epshead(i,j) = epshead(j,i)
2253 sig0head(i,j) = sig0head(j,i)
2256 wqdip(k,i,j) = wqdip(k,j,i)
2259 wquad(i,j) = wquad(j,i)
2260 epsintab(i,j) = epsintab(j,i)
2261 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2266 !--------------------- GBV potential -----------------------------------
2268 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2269 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2270 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2271 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2273 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2274 write (iout,'(a/)') 'The epsilon array:'
2275 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2276 write (iout,'(/a)') 'One-body parameters:'
2277 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2278 's||/s_|_^2',' chip ',' alph '
2279 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2280 sigii(i),chip(i),alp(i),i=1,ntyp)
2283 write(iout,*)"Wrong ipot"
2289 !-----------------------------------------------------------------------
2290 ! Calculate the "working" parameters of SC interactions.
2292 !el from module energy - COMMON.INTERACT-------
2293 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2294 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2295 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2296 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2297 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2298 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2300 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2306 if (scelemode.eq.0) then
2315 sc_aa_tube_par(:)=0.0d0
2316 sc_bb_tube_par(:)=0.0d0
2318 !--------------------------------
2323 epslip(i,j)=epslip(j,i)
2326 if (scelemode.eq.0) then
2329 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2330 sigma(j,i)=sigma(i,j)
2331 rs0(i,j)=dwa16*sigma(i,j)
2336 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2337 'Working parameters of the SC interactions:',&
2338 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2343 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2345 ! print *,"SIGMA ZLA?",sigma(i,j)
2353 sigeps=dsign(1.0D0,epsij)
2355 aa_aq(i,j)=epsij*rrij*rrij
2356 print *,"ADASKO",epsij,rrij,expon
2357 bb_aq(i,j)=-sigeps*epsij*rrij
2358 aa_aq(j,i)=aa_aq(i,j)
2359 bb_aq(j,i)=bb_aq(i,j)
2360 epsijlip=epslip(i,j)
2361 sigeps=dsign(1.0D0,epsijlip)
2362 epsijlip=dabs(epsijlip)
2363 aa_lip(i,j)=epsijlip*rrij*rrij
2364 bb_lip(i,j)=-sigeps*epsijlip*rrij
2365 aa_lip(j,i)=aa_lip(i,j)
2366 bb_lip(j,i)=bb_lip(i,j)
2367 !C write(iout,*) aa_lip
2368 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2369 sigt1sq=sigma0(i)**2
2370 sigt2sq=sigma0(j)**2
2373 ratsig1=sigt2sq/sigt1sq
2374 ratsig2=1.0D0/ratsig1
2375 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2376 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2377 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2381 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2382 sigmaii(i,j)=rsum_max
2383 sigmaii(j,i)=rsum_max
2385 ! sigmaii(i,j)=r0(i,j)
2386 ! sigmaii(j,i)=r0(i,j)
2388 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2389 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2390 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2391 augm(i,j)=epsij*r_augm**(2*expon)
2392 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2399 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2400 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2401 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2406 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2407 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2408 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2409 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2410 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2411 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2412 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2413 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2414 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2415 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2416 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2417 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2418 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2419 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2420 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2421 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2430 read (isidep_nucl,*) ipot_nucl
2431 ! print *,"TU?!",ipot_nucl
2432 if (ipot_nucl.eq.1) then
2433 do i=1,ntyp_molec(2)
2434 do j=i,ntyp_molec(2)
2435 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2436 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2440 do i=1,ntyp_molec(2)
2441 do j=i,ntyp_molec(2)
2442 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2443 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2444 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2448 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2449 do i=1,ntyp_molec(2)
2450 do j=i,ntyp_molec(2)
2451 rrij=sigma_nucl(i,j)
2455 epsij=4*eps_nucl(i,j)
2456 sigeps=dsign(1.0D0,epsij)
2458 aa_nucl(i,j)=epsij*rrij*rrij
2459 bb_nucl(i,j)=-sigeps*epsij*rrij
2460 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2461 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2462 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2463 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2464 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2465 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2468 aa_nucl(i,j)=aa_nucl(j,i)
2469 bb_nucl(i,j)=bb_nucl(j,i)
2470 ael3_nucl(i,j)=ael3_nucl(j,i)
2471 ael6_nucl(i,j)=ael6_nucl(j,i)
2472 ael63_nucl(i,j)=ael63_nucl(j,i)
2473 ael32_nucl(i,j)=ael32_nucl(j,i)
2474 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2475 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2476 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2477 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2478 eps_nucl(i,j)=eps_nucl(j,i)
2479 sigma_nucl(i,j)=sigma_nucl(j,i)
2480 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2484 write(iout,*) "tube param"
2485 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2486 ccavtubpep,dcavtubpep,tubetranenepep
2487 sigmapeptube=sigmapeptube**6
2488 sigeps=dsign(1.0D0,epspeptube)
2489 epspeptube=dabs(epspeptube)
2490 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2491 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2492 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2494 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2495 ccavtub(i),dcavtub(i),tubetranene(i)
2496 sigmasctube=sigmasctube**6
2497 sigeps=dsign(1.0D0,epssctube)
2498 epssctube=dabs(epssctube)
2499 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2500 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2501 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2503 !-----------------READING SC BASE POTENTIALS-----------------------------
2504 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2505 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2506 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2507 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2508 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2509 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2510 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2511 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2512 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2513 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2514 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2515 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2516 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2517 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2518 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2519 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2522 do i=1,ntyp_molec(1)
2523 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2524 write (*,*) "Im in ", i, " ", j
2525 read(isidep_scbase,*) &
2526 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2527 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2528 write(*,*) "eps",eps_scbase(i,j)
2529 read(isidep_scbase,*) &
2530 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2531 chis_scbase(i,j,1),chis_scbase(i,j,2)
2532 read(isidep_scbase,*) &
2533 dhead_scbasei(i,j), &
2534 dhead_scbasej(i,j), &
2535 rborn_scbasei(i,j),rborn_scbasej(i,j)
2536 read(isidep_scbase,*) &
2537 (wdipdip_scbase(k,i,j),k=1,3), &
2538 (wqdip_scbase(k,i,j),k=1,2)
2539 read(isidep_scbase,*) &
2540 alphapol_scbase(i,j), &
2541 epsintab_scbase(i,j)
2544 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2545 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2547 do i=1,ntyp_molec(1)
2548 do j=1,ntyp_molec(2)-1
2549 epsij=eps_scbase(i,j)
2550 rrij=sigma_scbase(i,j)
2555 sigeps=dsign(1.0D0,epsij)
2557 aa_scbase(i,j)=epsij*rrij*rrij
2558 bb_scbase(i,j)=-sigeps*epsij*rrij
2561 !-----------------READING PEP BASE POTENTIALS-------------------
2562 allocate(eps_pepbase(ntyp_molec(2)))
2563 allocate(sigma_pepbase(ntyp_molec(2)))
2564 allocate(chi_pepbase(ntyp_molec(2),2))
2565 allocate(chipp_pepbase(ntyp_molec(2),2))
2566 allocate(alphasur_pepbase(4,ntyp_molec(2)))
2567 allocate(sigmap1_pepbase(ntyp_molec(2)))
2568 allocate(sigmap2_pepbase(ntyp_molec(2)))
2569 allocate(chis_pepbase(ntyp_molec(2),2))
2570 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
2573 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2574 write (*,*) "Im in ", i, " ", j
2575 read(isidep_pepbase,*) &
2576 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
2577 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
2578 write(*,*) "eps",eps_pepbase(j)
2579 read(isidep_pepbase,*) &
2580 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
2581 chis_pepbase(j,1),chis_pepbase(j,2)
2582 read(isidep_pepbase,*) &
2583 (wdipdip_pepbase(k,j),k=1,3)
2585 allocate(aa_pepbase(ntyp_molec(2)))
2586 allocate(bb_pepbase(ntyp_molec(2)))
2588 do j=1,ntyp_molec(2)-1
2589 epsij=eps_pepbase(j)
2590 rrij=sigma_pepbase(j)
2595 sigeps=dsign(1.0D0,epsij)
2597 aa_pepbase(j)=epsij*rrij*rrij
2598 bb_pepbase(j)=-sigeps*epsij*rrij
2600 !--------------READING SC PHOSPHATE-------------------------------------
2601 allocate(eps_scpho(ntyp_molec(1)))
2602 allocate(sigma_scpho(ntyp_molec(1)))
2603 allocate(chi_scpho(ntyp_molec(1),2))
2604 allocate(chipp_scpho(ntyp_molec(1),2))
2605 allocate(alphasur_scpho(4,ntyp_molec(1)))
2606 allocate(sigmap1_scpho(ntyp_molec(1)))
2607 allocate(sigmap2_scpho(ntyp_molec(1)))
2608 allocate(chis_scpho(ntyp_molec(1),2))
2609 allocate(wqq_scpho(ntyp_molec(1)))
2610 allocate(wqdip_scpho(2,ntyp_molec(1)))
2611 allocate(alphapol_scpho(ntyp_molec(1)))
2612 allocate(epsintab_scpho(ntyp_molec(1)))
2613 allocate(dhead_scphoi(ntyp_molec(1)))
2614 allocate(rborn_scphoi(ntyp_molec(1)))
2615 allocate(rborn_scphoj(ntyp_molec(1)))
2616 allocate(alphi_scpho(ntyp_molec(1)))
2620 do j=1,ntyp_molec(1) ! without U then we will take T for U
2621 write (*,*) "Im in scpho ", i, " ", j
2622 read(isidep_scpho,*) &
2623 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
2624 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
2625 write(*,*) "eps",eps_scpho(j)
2626 read(isidep_scpho,*) &
2627 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
2628 chis_scpho(j,1),chis_scpho(j,2)
2629 read(isidep_scpho,*) &
2630 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
2631 read(isidep_scpho,*) &
2632 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
2636 allocate(aa_scpho(ntyp_molec(1)))
2637 allocate(bb_scpho(ntyp_molec(1)))
2639 do j=1,ntyp_molec(1)
2646 sigeps=dsign(1.0D0,epsij)
2648 aa_scpho(j)=epsij*rrij*rrij
2649 bb_scpho(j)=-sigeps*epsij*rrij
2653 read(isidep_peppho,*) &
2654 eps_peppho,sigma_peppho
2655 read(isidep_peppho,*) &
2656 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
2657 read(isidep_peppho,*) &
2658 (wqdip_peppho(k),k=1,2)
2666 sigeps=dsign(1.0D0,epsij)
2668 aa_peppho=epsij*rrij*rrij
2669 bb_peppho=-sigeps*epsij*rrij
2672 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
2677 ! Define the SC-p interaction constants (hard-coded; old style)
2680 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
2682 ! aad(i,1)=0.3D0*4.0D0**12
2683 ! Following line for constants currently implemented
2684 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
2685 aad(i,1)=1.5D0*4.0D0**12
2686 ! aad(i,1)=0.17D0*5.6D0**12
2688 ! "Soft" SC-p repulsion
2690 ! Following line for constants currently implemented
2691 ! aad(i,1)=0.3D0*4.0D0**6
2692 ! "Hard" SC-p repulsion
2693 bad(i,1)=3.0D0*4.0D0**6
2694 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
2703 ! 8/9/01 Read the SC-p interaction constants from file
2706 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
2709 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
2710 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
2711 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
2712 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
2716 write (iout,*) "Parameters of SC-p interactions:"
2718 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
2719 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
2724 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
2726 do i=1,ntyp_molec(2)
2727 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
2729 do i=1,ntyp_molec(2)
2730 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
2731 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
2733 r0pp=1.12246204830937298142*5.16158
2739 ! Define the constants of the disulfide bridge
2743 ! Old arbitrary potential - commented out.
2748 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
2749 ! energy surface of diethyl disulfide.
2750 ! A. Liwo and U. Kozlowska, 11/24/03
2767 write (iout,'(/a)') "Disulfide bridge parameters:"
2768 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
2769 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
2770 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
2771 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
2774 if (shield_mode.gt.0) then
2775 pi=4.0D0*datan(1.0D0)
2776 !C VSolvSphere the volume of solving sphere
2778 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
2779 !C there will be no distinction between proline peptide group and normal peptide
2780 !C group in case of shielding parameters
2781 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
2782 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
2783 write (iout,*) VSolvSphere,VSolvSphere_div
2784 !C long axis of side chain
2786 long_r_sidechain(i)=vbldsc0(1,i)
2787 ! if (scelemode.eq.0) then
2788 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
2789 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
2791 ! short_r_sidechain(i)=sigma(i,i)
2793 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
2800 111 write (iout,*) "Error reading bending energy parameters."
2802 112 write (iout,*) "Error reading rotamer energy parameters."
2804 113 write (iout,*) "Error reading torsional energy parameters."
2806 114 write (iout,*) "Error reading double torsional energy parameters."
2808 115 write (iout,*) &
2809 "Error reading cumulant (multibody energy) parameters."
2811 116 write (iout,*) "Error reading electrostatic energy parameters."
2813 117 write (iout,*) "Error reading side chain interaction parameters."
2815 118 write (iout,*) "Error reading SCp interaction parameters."
2817 119 write (iout,*) "Error reading SCCOR parameters"
2820 call MPI_Finalize(Ierror)
2824 end subroutine parmread
2826 !-----------------------------------------------------------------------------
2828 !-----------------------------------------------------------------------------
2829 subroutine printmat(ldim,m,n,iout,key,a)
2832 character(len=3),dimension(n) :: key
2833 real(kind=8),dimension(ldim,n) :: a
2835 integer :: i,j,k,m,iout,nlim
2839 write (iout,1000) (key(k),k=i,nlim)
2841 1000 format (/5x,8(6x,a3))
2842 1020 format (/80(1h-)/)
2844 write (iout,1010) key(j),(a(j,k),k=i,nlim)
2847 1010 format (a3,2x,8(f9.4))
2849 end subroutine printmat
2850 !-----------------------------------------------------------------------------
2852 !-----------------------------------------------------------------------------
2854 ! Read the PDB file and convert the peptide geometry into virtual-chain
2857 use energy_data, only: itype,istype
2861 ! use control, only: rescode,sugarcode
2862 ! implicit real*8 (a-h,o-z)
2863 ! include 'DIMENSIONS'
2864 ! include 'COMMON.LOCAL'
2865 ! include 'COMMON.VAR'
2866 ! include 'COMMON.CHAIN'
2867 ! include 'COMMON.INTERACT'
2868 ! include 'COMMON.IOUNITS'
2869 ! include 'COMMON.GEO'
2870 ! include 'COMMON.NAMES'
2871 ! include 'COMMON.CONTROL'
2872 ! include 'COMMON.DISTFIT'
2873 ! include 'COMMON.SETUP'
2874 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
2876 logical :: lprn=.true.,fail
2877 real(kind=8),dimension(3) :: e1,e2,e3
2878 real(kind=8) :: dcj,efree_temp
2879 character(len=3) :: seq,res
2880 character(len=5) :: atom
2881 character(len=80) :: card
2882 real(kind=8),dimension(3,20) :: sccor
2883 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
2884 integer :: isugar,molecprev,firstion
2885 character*1 :: sugar
2887 real(kind=8),dimension(3) :: ccc
2889 integer,dimension(2,maxres/3) :: hfrag_alloc
2890 integer,dimension(4,maxres/3) :: bfrag_alloc
2891 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
2892 real(kind=8),dimension(:,:), allocatable :: c_temporary
2893 integer,dimension(:,:) , allocatable :: itype_temporary
2894 integer,dimension(:),allocatable :: istype_temp
2901 ! write (2,*) "UNRES_PDB",unres_pdb
2921 !-----------------------------
2922 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
2923 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
2924 if(.not. allocated(istype)) allocate(istype(maxres))
2926 read (ipdbin,'(a80)',end=10) card
2927 write (iout,'(a)') card
2928 if (card(:5).eq.'HELIX') then
2931 read(card(22:25),*) hfrag(1,nhfrag)
2932 read(card(34:37),*) hfrag(2,nhfrag)
2934 if (card(:5).eq.'SHEET') then
2937 read(card(24:26),*) bfrag(1,nbfrag)
2938 read(card(35:37),*) bfrag(2,nbfrag)
2939 !rc----------------------------------------
2940 !rc to be corrected !!!
2941 bfrag(3,nbfrag)=bfrag(1,nbfrag)
2942 bfrag(4,nbfrag)=bfrag(2,nbfrag)
2943 !rc----------------------------------------
2945 if (card(:3).eq.'END') then
2947 else if (card(:3).eq.'TER') then
2952 itype(ires_old,molecule)=ntyp1_molec(molecule)
2953 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
2954 nres_molec(molecule)=nres_molec(molecule)+2
2956 ! write (iout,*) "Chain ended",ires,ishift,ires_old
2959 dc(j,ires)=sccor(j,iii)
2962 call sccenter(ires,iii,sccor)
2968 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
2969 ! Fish out the ATOM cards.
2970 ! write(iout,*) 'card',card(1:20)
2971 if (index(card(1:4),'ATOM').gt.0) then
2972 read (card(12:16),*) atom
2973 ! write (iout,*) "! ",atom," !",ires
2974 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
2975 read (card(23:26),*) ires
2976 read (card(18:20),'(a3)') res
2977 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
2978 ! & " ires_old",ires_old
2979 ! write (iout,*) "ishift",ishift," ishift1",ishift1
2980 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
2981 if (ires-ishift+ishift1.ne.ires_old) then
2982 ! Calculate the CM of the preceding residue.
2983 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
2985 ! write (iout,*) "Calculating sidechain center iii",iii
2988 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
2991 call sccenter(ires_old,iii,sccor)
2995 ! Start new residue.
2996 if (res.eq.'Cl-' .or. res.eq.'Na+') then
2999 else if (ibeg.eq.1) then
3000 write (iout,*) "BEG ires",ires
3002 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3005 nres_molec(molecule)=nres_molec(molecule)+1
3007 ires=ires-ishift+ishift1
3009 ! write (iout,*) "ishift",ishift," ires",ires,&
3010 ! " ires_old",ires_old
3012 else if (ibeg.eq.2) then
3014 ishift=-ires_old+ires-1 !!!!!
3015 ishift1=ishift1-1 !!!!!
3016 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3017 ires=ires-ishift+ishift1
3018 ! print *,ires,ishift,ishift1
3022 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3023 ires=ires-ishift+ishift1
3026 ! print *,'atom',ires,atom
3027 if (res.eq.'ACE' .or. res.eq.'NHE') then
3030 if (atom.eq.'CA '.or.atom.eq.'N ') then
3032 itype(ires,molecule)=rescode(ires,res,0,molecule)
3034 ! nres_molec(molecule)=nres_molec(molecule)+1
3037 itype(ires,molecule)=rescode(ires,res(2:3),0,molecule)
3038 ! nres_molec(molecule)=nres_molec(molecule)+1
3039 read (card(19:19),'(a1)') sugar
3040 isugar=sugarcode(sugar,ires)
3041 ! if (ibeg.eq.1) then
3045 ! print *,"ires=",ires,istype(ires)
3051 ires=ires-ishift+ishift1
3053 ! write (iout,*) "ires_old",ires_old," ires",ires
3054 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3057 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3058 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3059 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3060 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3061 ! print *,ires,ishift,ishift1
3062 ! write (iout,*) "backbone ",atom
3064 write (iout,'(2i3,2x,a,3f8.3)') &
3065 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3068 nres_molec(molecule)=nres_molec(molecule)+1
3070 sccor(j,iii)=c(j,ires)
3072 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3073 atom.eq."C2'" .or. atom.eq."C3'" &
3074 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3075 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3076 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3077 ! print *,ires,ishift,ishift1
3081 ! sccor(j,iii)=c(j,ires)
3084 c(j,ires)=c(j,ires)+ccc(j)/5.0
3086 print *,counter,molecule
3087 if (counter.eq.5) then
3089 nres_molec(molecule)=nres_molec(molecule)+1
3092 ! sccor(j,iii)=c(j,ires)
3096 ! print *, "ATOM",atom(1:3)
3097 else if (atom.eq."C5'") then
3098 read (card(19:19),'(a1)') sugar
3099 isugar=sugarcode(sugar,ires)
3104 ! print *,ires,istype(ires)
3108 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3109 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3110 nres_molec(molecule)=nres_molec(molecule)+1
3111 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3115 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3117 else if ((atom.eq."C1'").and.unres_pdb) then
3119 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3120 ! write (*,*) card(23:27),ires,itype(ires,1)
3121 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3122 atom.ne.'N' .and. atom.ne.'C' .and. &
3123 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3124 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3125 .and. atom.ne.'P '.and. &
3126 atom(1:1).ne.'H' .and. &
3127 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3129 ! write (iout,*) "sidechain ",atom
3130 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3131 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3132 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3134 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3137 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3138 if (firstion.eq.0) then
3142 dc(j,ires)=sccor(j,iii)
3145 call sccenter(ires,iii,sccor)
3148 read (card(12:16),*) atom
3149 print *,"HETATOM", atom
3150 read (card(18:20),'(a3)') res
3151 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3152 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3153 .or.(atom(1:2).eq.'K ')) &
3156 if (molecule.ne.5) molecprev=molecule
3158 nres_molec(molecule)=nres_molec(molecule)+1
3159 print *,"HERE",nres_molec(molecule)
3161 itype(ires,molecule)=rescode(ires,res,0,molecule)
3162 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3166 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3167 if (ires.eq.0) return
3168 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3171 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3173 nres_molec(molecule)=nres_molec(molecule)-2
3174 print *,'I have',nres, nres_molec(:)
3176 do k=1,4 ! ions are without dummy
3177 if (nres_molec(k).eq.0) cycle
3179 ! write (iout,*) i,itype(i,1)
3180 ! if (itype(i,1).eq.ntyp1) then
3181 ! write (iout,*) "dummy",i,itype(i,1)
3183 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3184 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3188 if (itype(i,k).eq.ntyp1_molec(k)) then
3189 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3190 if (itype(i+2,k).eq.0) then
3191 ! print *,"masz sieczke"
3193 if (itype(i+2,j).ne.0) then
3195 itype(i+1,j)=ntyp1_molec(j)
3196 nres_molec(k)=nres_molec(k)-1
3197 nres_molec(j)=nres_molec(j)+1
3203 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3204 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3205 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3206 ! if (unres_pdb) then
3207 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3208 ! print *,i,'tu dochodze'
3209 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3217 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3221 dcj=(c(j,i-2)-c(j,i-3))/2.0
3222 if (dcj.eq.0) dcj=1.23591524223
3227 else !itype(i+1,1).eq.ntyp1
3228 ! if (unres_pdb) then
3229 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3230 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3237 c(j,i)=c(j,i+1)-1.9d0*e2(j)
3241 dcj=(c(j,i+3)-c(j,i+2))/2.0
3242 if (dcj.eq.0) dcj=1.23591524223
3247 endif !itype(i+1,1).eq.ntyp1
3248 endif !itype.eq.ntyp1
3252 ! Calculate the CM of the last side chain.
3256 dc(j,ires)=sccor(j,iii)
3259 call sccenter(ires,iii,sccor)
3265 ! print *,"molecule",molecule
3266 if ((itype(nres,1).ne.10)) then
3268 if (molecule.eq.5) molecule=molecprev
3269 itype(nres,molecule)=ntyp1_molec(molecule)
3270 nres_molec(molecule)=nres_molec(molecule)+1
3271 ! if (unres_pdb) then
3272 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3273 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3280 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3284 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3285 c(j,nres)=c(j,nres-1)+dcj
3286 c(j,2*nres)=c(j,nres)
3290 ! print *,'I have',nres, nres_molec(:)
3292 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3294 if (nres.ne.nres0) then
3295 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3297 stop "Error nres value in WHAM input"
3300 !---------------------------------
3301 !el reallocate tables
3304 ! hfrag_alloc(j,i)=hfrag(j,i)
3307 ! bfrag_alloc(j,i)=bfrag(j,i)
3313 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3314 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3315 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3316 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3320 ! hfrag(j,i)=hfrag_alloc(j,i)
3325 ! bfrag(j,i)=bfrag_alloc(j,i)
3328 !el end reallocate tables
3329 !---------------------------------
3337 c(j,2*nres)=c(j,nres)
3340 if (itype(1,1).eq.ntyp1) then
3344 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3345 call refsys(2,3,4,e1,e2,e3,fail)
3352 c(j,1)=c(j,2)-1.9d0*e2(j)
3356 dcj=(c(j,4)-c(j,3))/2.0
3362 ! First lets assign correct dummy to correct type of chain
3364 if (itype(1,1).eq.ntyp1) then
3365 if (itype(2,1).eq.0) then
3367 if (itype(2,j).ne.0) then
3369 itype(1,j)=ntyp1_molec(j)
3370 nres_molec(1)=nres_molec(1)-1
3371 nres_molec(j)=nres_molec(j)+1
3378 print *,'I have',nres, nres_molec(:)
3380 ! Copy the coordinates to reference coordinates
3386 ! Calculate internal coordinates.
3388 write (iout,'(/a)') &
3389 "Cartesian coordinates of the reference structure"
3390 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3391 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3393 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3394 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3395 (c(j,ires+nres),j=1,3)
3398 ! znamy już nres wiec mozna alokowac tablice
3399 ! Calculate internal coordinates.
3400 if(me.eq.king.or..not.out1file)then
3401 write (iout,'(a)') &
3402 "Backbone and SC coordinates as read from the PDB"
3404 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
3405 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3406 (c(j,nres+ires),j=1,3)
3409 ! NOW LETS ROCK! SORTING
3410 allocate(c_temporary(3,2*nres))
3411 allocate(itype_temporary(nres,5))
3412 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3413 if (.not.allocated(istype)) write(iout,*) &
3414 "SOMETHING WRONG WITH ISTYTPE"
3415 allocate(istype_temp(nres))
3416 itype_temporary(:,:)=0
3420 if (itype(i,k).ne.0) then
3422 c_temporary(j,seqalingbegin)=c(j,i)
3423 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3426 itype_temporary(seqalingbegin,k)=itype(i,k)
3427 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3428 istype_temp(seqalingbegin)=istype(i)
3429 molnum(seqalingbegin)=k
3430 seqalingbegin=seqalingbegin+1
3436 c(j,i)=c_temporary(j,i)
3441 itype(i,k)=itype_temporary(i,k)
3442 istype(i)=istype_temp(i)
3445 ! if (itype(1,1).eq.ntyp1) then
3448 ! if (unres_pdb) then
3449 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3450 ! call refsys(2,3,4,e1,e2,e3,fail)
3457 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3461 ! dcj=(c(j,4)-c(j,3))/2.0
3463 ! c(j,nres+1)=c(j,1)
3469 write (iout,'(/a)') &
3470 "Cartesian coordinates of the reference structure after sorting"
3471 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3472 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3474 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3475 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3476 (c(j,ires+nres),j=1,3)
3480 ! print *,seqalingbegin,nres
3481 if(.not.allocated(vbld)) then
3482 allocate(vbld(2*nres))
3487 if(.not.allocated(vbld_inv)) then
3488 allocate(vbld_inv(2*nres))
3494 if(.not.allocated(theta)) then
3495 allocate(theta(nres+2))
3499 if(.not.allocated(phi)) allocate(phi(nres+2))
3500 if(.not.allocated(alph)) allocate(alph(nres+2))
3501 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3502 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3503 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3504 if(.not.allocated(costtab)) allocate(costtab(nres))
3505 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3506 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3507 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3508 if(.not.allocated(xxref)) allocate(xxref(nres))
3509 if(.not.allocated(yyref)) allocate(yyref(nres))
3510 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3511 if(.not.allocated(dc_norm)) then
3512 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3513 allocate(dc_norm(3,0:2*nres+2))
3517 call int_from_cart(.true.,.false.)
3518 call sc_loc_geom(.false.)
3520 thetaref(i)=theta(i)
3530 dc(j,i)=c(j,i+1)-c(j,i)
3531 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3536 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3537 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3539 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3543 ! Copy the coordinates to reference coordinates
3544 ! Splits to single chain if occurs
3548 ! cref(j,i,cou)=c(j,i)
3552 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
3553 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
3554 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
3555 !-----------------------------
3559 write (iout,*) "symetr", symetr
3562 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3564 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
3567 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
3572 cref(j,i,cou)=c(j,i)
3573 cref(j,i+nres,cou)=c(j,i+nres)
3575 chain_rep(j,lll,kkk)=c(j,i)
3576 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
3580 write (iout,*) chain_length
3581 if (chain_length.eq.0) chain_length=nres
3583 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
3584 chain_rep(j,chain_length+nres,symetr) &
3585 =chain_rep(j,chain_length+nres,1)
3588 ! write (iout,*) "spraw lancuchy",chain_length,symetr
3590 ! do kkk=1,chain_length
3591 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
3595 ! makes copy of chains
3596 write (iout,*) "symetr", symetr
3601 if (symetr.gt.1) then
3608 write(iout,*) (tabperm(i,kkk),kkk=1,4)
3614 write (iout,*) i,icha
3615 do lll=1,chain_length
3617 if (cou.le.nres) then
3619 kupa=mod(lll,chain_length)
3620 iprzes=(kkk-1)*chain_length+lll
3621 if (kupa.eq.0) kupa=chain_length
3622 write (iout,*) "kupa", kupa
3623 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
3624 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
3631 !-koniec robienia kopii
3634 write (iout,*) "nowa struktura", nperm
3636 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
3638 cref(3,i,kkk),cref(1,nres+i,kkk),&
3639 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
3641 100 format (//' alpha-carbon coordinates ',&
3642 ' centroid coordinates'/ &
3643 ' ', 6X,'X',11X,'Y',11X,'Z', &
3644 10X,'X',11X,'Y',11X,'Z')
3645 110 format (a,'(',i5,')',6f12.5)
3651 bfrag(i,j)=bfrag(i,j)-ishift
3657 hfrag(i,j)=hfrag(i,j)-ishift
3663 end subroutine readpdb
3664 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
3665 !-----------------------------------------------------------------------------
3667 !-----------------------------------------------------------------------------
3668 subroutine read_control
3682 use random, only: random_init
3683 ! implicit real*8 (a-h,o-z)
3684 ! include 'DIMENSIONS'
3686 use prng, only:prng_restart
3688 logical :: OKRandom!, prng_restart
3691 ! include 'COMMON.IOUNITS'
3692 ! include 'COMMON.TIME1'
3693 ! include 'COMMON.THREAD'
3694 ! include 'COMMON.SBRIDGE'
3695 ! include 'COMMON.CONTROL'
3696 ! include 'COMMON.MCM'
3697 ! include 'COMMON.MAP'
3698 ! include 'COMMON.HEADER'
3699 ! include 'COMMON.CSA'
3700 ! include 'COMMON.CHAIN'
3701 ! include 'COMMON.MUCA'
3702 ! include 'COMMON.MD'
3703 ! include 'COMMON.FFIELD'
3704 ! include 'COMMON.INTERACT'
3705 ! include 'COMMON.SETUP'
3706 !el integer :: KDIAG,ICORFL,IXDR
3707 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
3708 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
3709 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
3710 ! character(len=80) :: ucase
3711 character(len=640) :: controlcard
3713 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
3719 read (INP,'(a)') titel
3720 call card_concat(controlcard,.true.)
3721 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
3722 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
3723 call reada(controlcard,'SEED',seed,0.0D0)
3724 call random_init(seed)
3725 ! Set up the time limit (caution! The time must be input in minutes!)
3726 read_cart=index(controlcard,'READ_CART').gt.0
3727 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
3728 call readi(controlcard,'SYM',symetr,1)
3729 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
3730 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
3731 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
3732 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
3733 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
3734 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
3735 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
3736 call reada(controlcard,'DRMS',drms,0.1D0)
3737 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3738 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
3739 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
3740 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
3741 write (iout,'(a,f10.1)')'DRMS = ',drms
3742 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
3743 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
3745 call readi(controlcard,'NZ_START',nz_start,0)
3746 call readi(controlcard,'NZ_END',nz_end,0)
3747 call readi(controlcard,'IZ_SC',iz_sc,0)
3748 timlim=60.0D0*timlim
3749 safety = 60.0d0*safety
3752 call reada(controlcard,"T_BATH",t_bath,300.0d0)
3753 !C SHIELD keyword sets if the shielding effect of side-chains is used
3754 !C 0 denots no shielding is used all peptide are equally despite the
3755 !C solvent accesible area
3756 !C 1 the newly introduced function
3757 !C 2 reseved for further possible developement
3758 call readi(controlcard,'SHIELD',shield_mode,0)
3759 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3760 write(iout,*) "shield_mode",shield_mode
3761 !C Varibles set size of box
3762 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
3763 protein=index(controlcard,"PROTEIN").gt.0
3764 ions=index(controlcard,"IONS").gt.0
3765 nucleic=index(controlcard,"NUCLEIC").gt.0
3766 write (iout,*) "with_theta_constr ",with_theta_constr
3767 AFMlog=(index(controlcard,'AFM'))
3768 selfguide=(index(controlcard,'SELFGUIDE'))
3769 print *,'AFMlog',AFMlog,selfguide,"KUPA"
3770 call readi(controlcard,'GENCONSTR',genconstr,0)
3771 call reada(controlcard,'BOXX',boxxsize,100.0d0)
3772 call reada(controlcard,'BOXY',boxysize,100.0d0)
3773 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
3774 call readi(controlcard,'TUBEMOD',tubemode,0)
3775 print *,"SCELE",scelemode
3776 call readi(controlcard,"SCELEMODE",scelemode,0)
3777 print *,"SCELE",scelemode
3779 ! elemode = 0 is orignal UNRES electrostatics
3780 ! elemode = 1 is "Momo" potentials in progress
3781 ! elemode = 2 is in development EVALD
3782 write (iout,*) TUBEmode,"TUBEMODE"
3783 if (TUBEmode.gt.0) then
3784 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
3785 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
3786 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
3787 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
3788 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
3789 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
3790 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
3791 buftubebot=bordtubebot+tubebufthick
3792 buftubetop=bordtubetop-tubebufthick
3795 ! CUTOFFF ON ELECTROSTATICS
3796 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
3797 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
3798 write(iout,*) "R_CUT_ELE=",r_cut_ele
3799 ! Lipidic parameters
3800 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
3801 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
3802 if (lipthick.gt.0.0d0) then
3803 bordliptop=(boxzsize+lipthick)/2.0
3804 bordlipbot=bordliptop-lipthick
3805 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
3806 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
3807 buflipbot=bordlipbot+lipbufthick
3808 bufliptop=bordliptop-lipbufthick
3809 if ((lipbufthick*2.0d0).gt.lipthick) &
3810 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
3811 endif !lipthick.gt.0
3812 write(iout,*) "bordliptop=",bordliptop
3813 write(iout,*) "bordlipbot=",bordlipbot
3814 write(iout,*) "bufliptop=",bufliptop
3815 write(iout,*) "buflipbot=",buflipbot
3816 write (iout,*) "SHIELD MODE",shield_mode
3818 !C-------------------------
3819 minim=(index(controlcard,'MINIMIZE').gt.0)
3820 dccart=(index(controlcard,'CART').gt.0)
3821 overlapsc=(index(controlcard,'OVERLAP').gt.0)
3822 overlapsc=.not.overlapsc
3823 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
3824 searchsc=.not.searchsc
3825 sideadd=(index(controlcard,'SIDEADD').gt.0)
3826 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
3827 outpdb=(index(controlcard,'PDBOUT').gt.0)
3828 outmol2=(index(controlcard,'MOL2OUT').gt.0)
3829 pdbref=(index(controlcard,'PDBREF').gt.0)
3830 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
3831 indpdb=index(controlcard,'PDBSTART')
3832 extconf=(index(controlcard,'EXTCONF').gt.0)
3833 call readi(controlcard,'IPRINT',iprint,0)
3834 call readi(controlcard,'MAXGEN',maxgen,10000)
3835 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
3836 call readi(controlcard,"KDIAG",kdiag,0)
3837 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
3838 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
3839 write (iout,*) "RESCALE_MODE",rescale_mode
3840 split_ene=index(controlcard,'SPLIT_ENE').gt.0
3841 if (index(controlcard,'REGULAR').gt.0.0D0) then
3842 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3846 if (index(controlcard,'CHECKGRAD').gt.0) then
3848 if (index(controlcard,'CART').gt.0) then
3850 elseif (index(controlcard,'CARINT').gt.0) then
3855 elseif (index(controlcard,'THREAD').gt.0) then
3857 call readi(controlcard,'THREAD',nthread,0)
3858 if (nthread.gt.0) then
3859 call reada(controlcard,'WEIDIS',weidis,0.1D0)
3862 write (iout,'(a)')'A number has to follow the THREAD keyword.'
3863 stop 'Error termination in Read_Control.'
3865 else if (index(controlcard,'MCMA').gt.0) then
3867 else if (index(controlcard,'MCEE').gt.0) then
3869 else if (index(controlcard,'MULTCONF').gt.0) then
3871 else if (index(controlcard,'MAP').gt.0) then
3873 call readi(controlcard,'MAP',nmap,0)
3874 else if (index(controlcard,'CSA').gt.0) then
3876 !rc else if (index(controlcard,'ZSCORE').gt.0) then
3878 !rc ZSCORE is rm from UNRES, modecalc=9 is available
3881 !fcm else if (index(controlcard,'MCMF').gt.0) then
3883 else if (index(controlcard,'SOFTREG').gt.0) then
3885 else if (index(controlcard,'CHECK_BOND').gt.0) then
3887 else if (index(controlcard,'TEST').gt.0) then
3889 else if (index(controlcard,'MD').gt.0) then
3891 else if (index(controlcard,'RE ').gt.0) then
3895 lmuca=index(controlcard,'MUCA').gt.0
3896 call readi(controlcard,'MUCADYN',mucadyn,0)
3897 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
3898 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
3900 write (iout,*) 'MUCADYN=',mucadyn
3901 write (iout,*) 'MUCASMOOTH=',muca_smooth
3904 iscode=index(controlcard,'ONE_LETTER')
3905 indphi=index(controlcard,'PHI')
3906 indback=index(controlcard,'BACK')
3907 iranconf=index(controlcard,'RAND_CONF')
3908 i2ndstr=index(controlcard,'USE_SEC_PRED')
3909 gradout=index(controlcard,'GRADOUT').gt.0
3910 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
3911 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
3912 if (me.eq.king .or. .not.out1file ) &
3913 write (iout,*) "DISTCHAINMAX",distchainmax
3915 if(me.eq.king.or..not.out1file) &
3916 write (iout,'(2a)') diagmeth(kdiag),&
3917 ' routine used to diagonalize matrices.'
3918 if (shield_mode.gt.0) then
3919 pi=4.0D0*datan(1.0D0)
3920 !C VSolvSphere the volume of solving sphere
3922 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3923 !C there will be no distinction between proline peptide group and normal peptide
3924 !C group in case of shielding parameters
3925 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3926 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3927 write (iout,*) VSolvSphere,VSolvSphere_div
3928 !C long axis of side chain
3930 ! long_r_sidechain(i)=vbldsc0(1,i)
3931 ! short_r_sidechain(i)=sigma0(i)
3932 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3938 end subroutine read_control
3939 !-----------------------------------------------------------------------------
3940 subroutine read_REMDpar
3942 ! Read REMD settings
3949 use control_data, only:out1file
3950 ! implicit real*8 (a-h,o-z)
3951 ! include 'DIMENSIONS'
3952 ! include 'COMMON.IOUNITS'
3953 ! include 'COMMON.TIME1'
3954 ! include 'COMMON.MD'
3957 !el include 'COMMON.LANGEVIN'
3959 !el include 'COMMON.LANGEVIN.lang0'
3961 ! include 'COMMON.INTERACT'
3962 ! include 'COMMON.NAMES'
3963 ! include 'COMMON.GEO'
3964 ! include 'COMMON.REMD'
3965 ! include 'COMMON.CONTROL'
3966 ! include 'COMMON.SETUP'
3967 ! character(len=80) :: ucase
3968 character(len=320) :: controlcard
3969 character(len=3200) :: controlcard1
3970 integer :: iremd_m_total
3973 ! real(kind=8) :: var,ene
3975 if(me.eq.king.or..not.out1file) &
3976 write (iout,*) "REMD setup"
3978 call card_concat(controlcard,.true.)
3979 call readi(controlcard,"NREP",nrep,3)
3980 call readi(controlcard,"NSTEX",nstex,1000)
3981 call reada(controlcard,"RETMIN",retmin,10.0d0)
3982 call reada(controlcard,"RETMAX",retmax,1000.0d0)
3983 mremdsync=(index(controlcard,'SYNC').gt.0)
3984 call readi(controlcard,"NSYN",i_sync_step,100)
3985 restart1file=(index(controlcard,'REST1FILE').gt.0)
3986 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
3987 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
3988 if(max_cache_traj_use.gt.max_cache_traj) &
3989 max_cache_traj_use=max_cache_traj
3990 if(me.eq.king.or..not.out1file) then
3991 !d if (traj1file) then
3992 !rc caching is in testing - NTWX is not ignored
3993 !d write (iout,*) "NTWX value is ignored"
3994 !d write (iout,*) " trajectory is stored to one file by master"
3995 !d write (iout,*) " before exchange at NSTEX intervals"
3997 write (iout,*) "NREP= ",nrep
3998 write (iout,*) "NSTEX= ",nstex
3999 write (iout,*) "SYNC= ",mremdsync
4000 write (iout,*) "NSYN= ",i_sync_step
4001 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4004 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4005 if (index(controlcard,'TLIST').gt.0) then
4007 call card_concat(controlcard1,.true.)
4008 read(controlcard1,*) (remd_t(i),i=1,nrep)
4009 if(me.eq.king.or..not.out1file) &
4010 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4013 if (index(controlcard,'MLIST').gt.0) then
4015 call card_concat(controlcard1,.true.)
4016 read(controlcard1,*) (remd_m(i),i=1,nrep)
4017 if(me.eq.king.or..not.out1file) then
4018 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4021 iremd_m_total=iremd_m_total+remd_m(i)
4023 write (iout,*) 'Total number of replicas ',iremd_m_total
4026 if(me.eq.king.or..not.out1file) &
4027 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4029 end subroutine read_REMDpar
4030 !-----------------------------------------------------------------------------
4031 subroutine read_MDpar
4035 use control_data, only: r_cut,rlamb,out1file
4037 use geometry_data, only: pi
4039 ! implicit real*8 (a-h,o-z)
4040 ! include 'DIMENSIONS'
4041 ! include 'COMMON.IOUNITS'
4042 ! include 'COMMON.TIME1'
4043 ! include 'COMMON.MD'
4046 !el include 'COMMON.LANGEVIN'
4048 !el include 'COMMON.LANGEVIN.lang0'
4050 ! include 'COMMON.INTERACT'
4051 ! include 'COMMON.NAMES'
4052 ! include 'COMMON.GEO'
4053 ! include 'COMMON.SETUP'
4054 ! include 'COMMON.CONTROL'
4055 ! include 'COMMON.SPLITELE'
4056 ! character(len=80) :: ucase
4057 character(len=320) :: controlcard
4062 call card_concat(controlcard,.true.)
4063 call readi(controlcard,"NSTEP",n_timestep,1000000)
4064 call readi(controlcard,"NTWE",ntwe,100)
4065 call readi(controlcard,"NTWX",ntwx,1000)
4066 call reada(controlcard,"DT",d_time,1.0d-1)
4067 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4068 call reada(controlcard,"DAMAX",damax,1.0d1)
4069 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4070 call readi(controlcard,"LANG",lang,0)
4071 RESPA = index(controlcard,"RESPA") .gt. 0
4072 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4073 ntime_split0=ntime_split
4074 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4075 ntime_split0=ntime_split
4076 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4077 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4078 rest = index(controlcard,"REST").gt.0
4079 tbf = index(controlcard,"TBF").gt.0
4080 usampl = index(controlcard,"USAMPL").gt.0
4081 mdpdb = index(controlcard,"MDPDB").gt.0
4082 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4083 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4084 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4085 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4086 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4087 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4088 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4089 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4090 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4091 large = index(controlcard,"LARGE").gt.0
4092 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4093 rattle = index(controlcard,"RATTLE").gt.0
4094 preminim=(index(controlcard,'PREMINIM').gt.0)
4095 write (iout,*) "PREMINIM ",preminim
4096 dccart=(index(controlcard,'CART').gt.0)
4097 if (preminim) call read_minim
4098 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4104 if(me.eq.king.or..not.out1file) then
4106 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4108 write (iout,'(a)') "The units are:"
4109 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4110 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4111 " acceleration: angstrom/(48.9 fs)**2"
4112 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4114 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4115 write (iout,'(a60,f10.5,a)') &
4116 "Initial time step of numerical integration:",d_time,&
4118 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4120 write (iout,'(2a,i4,a)') &
4121 "A-MTS algorithm used; initial time step for fast-varying",&
4122 " short-range forces split into",ntime_split," steps."
4123 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4124 r_cut," lambda",rlamb
4126 write (iout,'(2a,f10.5)') &
4127 "Maximum acceleration threshold to reduce the time step",&
4128 "/increase split number:",damax
4129 write (iout,'(2a,f10.5)') &
4130 "Maximum predicted energy drift to reduce the timestep",&
4131 "/increase split number:",edriftmax
4132 write (iout,'(a60,f10.5)') &
4133 "Maximum velocity threshold to reduce velocities:",dvmax
4134 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4135 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4136 if (rattle) write (iout,'(a60)') &
4137 "Rattle algorithm used to constrain the virtual bonds"
4141 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4142 call reada(controlcard,"RWAT",rwat,1.4d0)
4143 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4144 surfarea=index(controlcard,"SURFAREA").gt.0
4145 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4146 if(me.eq.king.or..not.out1file)then
4147 write (iout,'(/a,$)') "Langevin dynamics calculation"
4149 write (iout,'(a/)') &
4150 " with direct integration of Langevin equations"
4151 else if (lang.eq.2) then
4152 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4153 else if (lang.eq.3) then
4154 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4155 else if (lang.eq.4) then
4156 write (iout,'(a/)') " in overdamped mode"
4158 write (iout,'(//a,i5)') &
4159 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4162 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4163 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4164 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4165 write (iout,'(a60,f10.5)') &
4166 "Scaling factor of the friction forces:",scal_fric
4167 if (surfarea) write (iout,'(2a,i10,a)') &
4168 "Friction coefficients will be scaled by solvent-accessible",&
4169 " surface area every",reset_fricmat," steps."
4171 ! Calculate friction coefficients and bounds of stochastic forces
4172 eta=6*pi*cPoise*etawat
4173 if(me.eq.king.or..not.out1file) &
4174 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4177 do j=1,5 !types of molecules
4178 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4179 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4181 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4182 do j=1,5 !types of molecules
4184 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4185 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4189 if(me.eq.king.or..not.out1file)then
4190 write (iout,'(/2a/)') &
4191 "Radii of site types and friction coefficients and std's of",&
4192 " stochastic forces of fully exposed sites"
4193 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4195 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4196 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4200 if(me.eq.king.or..not.out1file)then
4201 write (iout,'(a)') "Berendsen bath calculation"
4202 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4203 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4205 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4206 count_reset_moment," steps"
4208 write (iout,'(a,i10,a)') &
4209 "Velocities will be reset at random every",count_reset_vel,&
4213 if(me.eq.king.or..not.out1file) &
4214 write (iout,'(a31)') "Microcanonical mode calculation"
4216 if(me.eq.king.or..not.out1file)then
4217 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4219 write(iout,*) "MD running with constraints."
4220 write(iout,*) "Equilibration time ", eq_time, " mtus."
4221 write(iout,*) "Constraining ", nfrag," fragments."
4222 write(iout,*) "Length of each fragment, weight and q0:"
4224 write (iout,*) "Set of restraints #",iset
4226 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4227 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4229 write(iout,*) "constraints between ", npair, "fragments."
4230 write(iout,*) "constraint pairs, weights and q0:"
4232 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4233 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4235 write(iout,*) "angle constraints within ", nfrag_back,&
4236 "backbone fragments."
4237 write(iout,*) "fragment, weights:"
4239 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4240 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4241 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4244 iset=mod(kolor,nset)+1
4247 if(me.eq.king.or..not.out1file) &
4248 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4250 end subroutine read_MDpar
4251 !-----------------------------------------------------------------------------
4255 ! implicit real*8 (a-h,o-z)
4256 ! include 'DIMENSIONS'
4257 ! include 'COMMON.MAP'
4258 ! include 'COMMON.IOUNITS'
4259 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4260 character(len=80) :: mapcard !,ucase
4263 ! real(kind=8) :: var,ene
4266 read (inp,'(a)') mapcard
4267 mapcard=ucase(mapcard)
4268 if (index(mapcard,'PHI').gt.0) then
4270 else if (index(mapcard,'THE').gt.0) then
4272 else if (index(mapcard,'ALP').gt.0) then
4274 else if (index(mapcard,'OME').gt.0) then
4277 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4278 stop 'Error - illegal variable spec in MAP card.'
4280 call readi (mapcard,'RES1',res1(imap),0)
4281 call readi (mapcard,'RES2',res2(imap),0)
4282 if (res1(imap).eq.0) then
4283 res1(imap)=res2(imap)
4284 else if (res2(imap).eq.0) then
4285 res2(imap)=res1(imap)
4287 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4288 write (iout,'(a)') &
4289 'Error - illegal definition of variable group in MAP.'
4290 stop 'Error - illegal definition of variable group in MAP.'
4292 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4293 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4294 call readi(mapcard,'NSTEP',nstep(imap),0)
4295 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4296 write (iout,'(a)') &
4297 'Illegal boundary and/or step size specification in MAP.'
4298 stop 'Illegal boundary and/or step size specification in MAP.'
4302 end subroutine map_read
4303 !-----------------------------------------------------------------------------
4306 use control_data, only: vdisulf
4308 ! implicit real*8 (a-h,o-z)
4309 ! include 'DIMENSIONS'
4310 ! include 'COMMON.IOUNITS'
4311 ! include 'COMMON.GEO'
4312 ! include 'COMMON.CSA'
4313 ! include 'COMMON.BANK'
4314 ! include 'COMMON.CONTROL'
4315 ! character(len=80) :: ucase
4316 character(len=620) :: mcmcard
4318 ! integer :: ntf,ik,iw_pdb
4319 ! real(kind=8) :: var,ene
4321 call card_concat(mcmcard,.true.)
4323 call readi(mcmcard,'NCONF',nconf,50)
4324 call readi(mcmcard,'NADD',nadd,0)
4325 call readi(mcmcard,'JSTART',jstart,1)
4326 call readi(mcmcard,'JEND',jend,1)
4327 call readi(mcmcard,'NSTMAX',nstmax,500000)
4328 call readi(mcmcard,'N0',n0,1)
4329 call readi(mcmcard,'N1',n1,6)
4330 call readi(mcmcard,'N2',n2,4)
4331 call readi(mcmcard,'N3',n3,0)
4332 call readi(mcmcard,'N4',n4,0)
4333 call readi(mcmcard,'N5',n5,0)
4334 call readi(mcmcard,'N6',n6,10)
4335 call readi(mcmcard,'N7',n7,0)
4336 call readi(mcmcard,'N8',n8,0)
4337 call readi(mcmcard,'N9',n9,0)
4338 call readi(mcmcard,'N14',n14,0)
4339 call readi(mcmcard,'N15',n15,0)
4340 call readi(mcmcard,'N16',n16,0)
4341 call readi(mcmcard,'N17',n17,0)
4342 call readi(mcmcard,'N18',n18,0)
4344 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4346 call readi(mcmcard,'NDIFF',ndiff,2)
4347 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4348 call readi(mcmcard,'IS1',is1,1)
4349 call readi(mcmcard,'IS2',is2,8)
4350 call readi(mcmcard,'NRAN0',nran0,4)
4351 call readi(mcmcard,'NRAN1',nran1,2)
4352 call readi(mcmcard,'IRR',irr,1)
4353 call readi(mcmcard,'NSEED',nseed,20)
4354 call readi(mcmcard,'NTOTAL',ntotal,10000)
4355 call reada(mcmcard,'CUT1',cut1,2.0d0)
4356 call reada(mcmcard,'CUT2',cut2,5.0d0)
4357 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4358 call readi(mcmcard,'ICMAX',icmax,3)
4359 call readi(mcmcard,'IRESTART',irestart,0)
4360 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4363 call reada(mcmcard,'DELE',dele,20.0d0)
4364 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4365 call readi(mcmcard,'IREF',iref,0)
4366 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4367 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4368 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4369 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4370 write (iout,*) "NCONF_IN",nconf_in
4372 end subroutine csaread
4373 !-----------------------------------------------------------------------------
4377 use control_data, only: MaxMoveType
4380 ! implicit real*8 (a-h,o-z)
4381 ! include 'DIMENSIONS'
4382 ! include 'COMMON.MCM'
4383 ! include 'COMMON.MCE'
4384 ! include 'COMMON.IOUNITS'
4385 ! character(len=80) :: ucase
4386 character(len=320) :: mcmcard
4389 ! real(kind=8) :: var,ene
4391 call card_concat(mcmcard,.true.)
4392 call readi(mcmcard,'MAXACC',maxacc,100)
4393 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4394 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4395 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4396 call readi(mcmcard,'MAXREPM',maxrepm,200)
4397 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4398 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4399 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4400 call reada(mcmcard,'E_UP',e_up,5.0D0)
4401 call reada(mcmcard,'DELTE',delte,0.1D0)
4402 call readi(mcmcard,'NSWEEP',nsweep,5)
4403 call readi(mcmcard,'NSTEPH',nsteph,0)
4404 call readi(mcmcard,'NSTEPC',nstepc,0)
4405 call reada(mcmcard,'TMIN',tmin,298.0D0)
4406 call reada(mcmcard,'TMAX',tmax,298.0D0)
4407 call readi(mcmcard,'NWINDOW',nwindow,0)
4408 call readi(mcmcard,'PRINT_MC',print_mc,0)
4409 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4410 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4411 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4412 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4413 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4414 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4415 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4416 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4417 if (nwindow.gt.0) then
4418 allocate(winstart(nwindow)) !!el (maxres)
4419 allocate(winend(nwindow)) !!el
4420 allocate(winlen(nwindow)) !!el
4421 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4423 winlen(i)=winend(i)-winstart(i)+1
4426 if (tmax.lt.tmin) tmax=tmin
4427 if (tmax.eq.tmin) then
4431 if (nstepc.gt.0 .and. nsteph.gt.0) then
4432 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4433 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4435 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4436 ! Probabilities of different move types
4437 sumpro_type(0)=0.0D0
4438 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4439 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4440 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4441 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4442 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4443 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4444 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4446 print *,'i',i,' sumprotype',sumpro_type(i)
4447 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4448 print *,'i',i,' sumprotype',sumpro_type(i)
4451 end subroutine mcmread
4452 !-----------------------------------------------------------------------------
4453 subroutine read_minim
4456 ! implicit real*8 (a-h,o-z)
4457 ! include 'DIMENSIONS'
4458 ! include 'COMMON.MINIM'
4459 ! include 'COMMON.IOUNITS'
4460 ! character(len=80) :: ucase
4461 character(len=320) :: minimcard
4463 ! integer :: ntf,ik,iw_pdb
4464 ! real(kind=8) :: var,ene
4466 call card_concat(minimcard,.true.)
4467 call readi(minimcard,'MAXMIN',maxmin,2000)
4468 call readi(minimcard,'MAXFUN',maxfun,5000)
4469 call readi(minimcard,'MINMIN',minmin,maxmin)
4470 call readi(minimcard,'MINFUN',minfun,maxmin)
4471 call reada(minimcard,'TOLF',tolf,1.0D-2)
4472 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4473 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4474 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4475 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4476 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4477 'Options in energy minimization:'
4478 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4479 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4480 'MinMin:',MinMin,' MinFun:',MinFun,&
4481 ' TolF:',TolF,' RTolF:',RTolF
4483 end subroutine read_minim
4484 !-----------------------------------------------------------------------------
4485 subroutine openunits
4487 use MD_data, only: usampl
4490 use control_data, only:out1file
4491 use control, only: getenv_loc
4492 ! implicit real*8 (a-h,o-z)
4493 ! include 'DIMENSIONS'
4496 character(len=16) :: form,nodename
4497 integer :: nodelen,ierror,npos
4499 ! include 'COMMON.SETUP'
4500 ! include 'COMMON.IOUNITS'
4501 ! include 'COMMON.MD'
4502 ! include 'COMMON.CONTROL'
4503 integer :: lenpre,lenpot,lentmp !,ilen
4505 character(len=3) :: out1file_text !,ucase
4506 character(len=3) :: ll
4509 ! integer :: ntf,ik,iw_pdb
4510 ! real(kind=8) :: var,ene
4512 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4513 call getenv_loc("PREFIX",prefix)
4515 call getenv_loc("POT",pot)
4516 call getenv_loc("DIRTMP",tmpdir)
4517 call getenv_loc("CURDIR",curdir)
4518 call getenv_loc("OUT1FILE",out1file_text)
4519 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4520 out1file_text=ucase(out1file_text)
4521 if (out1file_text(1:1).eq."Y") then
4524 out1file=fg_rank.gt.0
4529 if (lentmp.gt.0) then
4530 write (*,'(80(1h!))')
4531 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4532 write (*,'(80(1h!))')
4533 write (*,*)"All output files will be on node /tmp directory."
4535 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4536 if (me.eq.king) then
4537 write (*,*) "The master node is ",nodename
4538 else if (fg_rank.eq.0) then
4539 write (*,*) "I am the CG slave node ",nodename
4541 write (*,*) "I am the FG slave node ",nodename
4544 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
4545 lenpre = lentmp+lenpre+1
4547 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
4548 ! Get the names and open the input files
4549 #if defined(WINIFL) || defined(WINPGI)
4550 open(1,file=pref_orig(:ilen(pref_orig))// &
4551 '.inp',status='old',readonly,shared)
4552 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4553 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4554 ! Get parameter filenames and open the parameter files.
4555 call getenv_loc('BONDPAR',bondname)
4556 open (ibond,file=bondname,status='old',readonly,shared)
4557 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4558 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
4559 call getenv_loc('THETPAR',thetname)
4560 open (ithep,file=thetname,status='old',readonly,shared)
4561 call getenv_loc('ROTPAR',rotname)
4562 open (irotam,file=rotname,status='old',readonly,shared)
4563 call getenv_loc('TORPAR',torname)
4564 open (itorp,file=torname,status='old',readonly,shared)
4565 call getenv_loc('TORDPAR',tordname)
4566 open (itordp,file=tordname,status='old',readonly,shared)
4567 call getenv_loc('FOURIER',fouriername)
4568 open (ifourier,file=fouriername,status='old',readonly,shared)
4569 call getenv_loc('ELEPAR',elename)
4570 open (ielep,file=elename,status='old',readonly,shared)
4571 call getenv_loc('SIDEPAR',sidename)
4572 open (isidep,file=sidename,status='old',readonly,shared)
4574 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4575 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
4576 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4577 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
4578 call getenv_loc('TORPAR_NUCL',torname_nucl)
4579 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
4580 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4581 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
4582 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4583 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
4586 #elif (defined CRAY) || (defined AIX)
4587 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4589 ! print *,"Processor",myrank," opened file 1"
4590 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4591 ! print *,"Processor",myrank," opened file 9"
4592 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4593 ! Get parameter filenames and open the parameter files.
4594 call getenv_loc('BONDPAR',bondname)
4595 open (ibond,file=bondname,status='old',action='read')
4596 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4597 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4599 ! print *,"Processor",myrank," opened file IBOND"
4600 call getenv_loc('THETPAR',thetname)
4601 open (ithep,file=thetname,status='old',action='read')
4602 ! print *,"Processor",myrank," opened file ITHEP"
4603 call getenv_loc('ROTPAR',rotname)
4604 open (irotam,file=rotname,status='old',action='read')
4605 ! print *,"Processor",myrank," opened file IROTAM"
4606 call getenv_loc('TORPAR',torname)
4607 open (itorp,file=torname,status='old',action='read')
4608 ! print *,"Processor",myrank," opened file ITORP"
4609 call getenv_loc('TORDPAR',tordname)
4610 open (itordp,file=tordname,status='old',action='read')
4611 ! print *,"Processor",myrank," opened file ITORDP"
4612 call getenv_loc('SCCORPAR',sccorname)
4613 open (isccor,file=sccorname,status='old',action='read')
4614 ! print *,"Processor",myrank," opened file ISCCOR"
4615 call getenv_loc('FOURIER',fouriername)
4616 open (ifourier,file=fouriername,status='old',action='read')
4617 ! print *,"Processor",myrank," opened file IFOURIER"
4618 call getenv_loc('ELEPAR',elename)
4619 open (ielep,file=elename,status='old',action='read')
4620 ! print *,"Processor",myrank," opened file IELEP"
4621 call getenv_loc('SIDEPAR',sidename)
4622 open (isidep,file=sidename,status='old',action='read')
4624 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4625 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4626 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4627 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4628 call getenv_loc('TORPAR_NUCL',torname_nucl)
4629 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4630 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4631 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4632 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4633 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4635 call getenv_loc('LIPTRANPAR',liptranname)
4636 open (iliptranpar,file=liptranname,status='old',action='read')
4637 call getenv_loc('TUBEPAR',tubename)
4638 open (itube,file=tubename,status='old',action='read')
4639 call getenv_loc('IONPAR',ionname)
4640 open (iion,file=ionname,status='old',action='read')
4642 ! print *,"Processor",myrank," opened file ISIDEP"
4643 ! print *,"Processor",myrank," opened parameter files"
4645 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
4646 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4647 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4648 ! Get parameter filenames and open the parameter files.
4649 call getenv_loc('BONDPAR',bondname)
4650 open (ibond,file=bondname,status='old')
4651 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4652 open (ibond_nucl,file=bondname_nucl,status='old')
4654 call getenv_loc('THETPAR',thetname)
4655 open (ithep,file=thetname,status='old')
4656 call getenv_loc('ROTPAR',rotname)
4657 open (irotam,file=rotname,status='old')
4658 call getenv_loc('TORPAR',torname)
4659 open (itorp,file=torname,status='old')
4660 call getenv_loc('TORDPAR',tordname)
4661 open (itordp,file=tordname,status='old')
4662 call getenv_loc('SCCORPAR',sccorname)
4663 open (isccor,file=sccorname,status='old')
4664 call getenv_loc('FOURIER',fouriername)
4665 open (ifourier,file=fouriername,status='old')
4666 call getenv_loc('ELEPAR',elename)
4667 open (ielep,file=elename,status='old')
4668 call getenv_loc('SIDEPAR',sidename)
4669 open (isidep,file=sidename,status='old')
4671 open (ithep_nucl,file=thetname_nucl,status='old')
4672 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4673 open (irotam_nucl,file=rotname_nucl,status='old')
4674 call getenv_loc('TORPAR_NUCL',torname_nucl)
4675 open (itorp_nucl,file=torname_nucl,status='old')
4676 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4677 open (itordp_nucl,file=tordname_nucl,status='old')
4678 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4679 open (isidep_nucl,file=sidename_nucl,status='old')
4681 call getenv_loc('LIPTRANPAR',liptranname)
4682 open (iliptranpar,file=liptranname,status='old')
4683 call getenv_loc('TUBEPAR',tubename)
4684 open (itube,file=tubename,status='old')
4685 call getenv_loc('IONPAR',ionname)
4686 open (iion,file=ionname,status='old')
4688 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
4690 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
4691 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
4692 ! Get parameter filenames and open the parameter files.
4693 call getenv_loc('BONDPAR',bondname)
4694 open (ibond,file=bondname,status='old',action='read')
4695 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
4696 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
4697 call getenv_loc('THETPAR',thetname)
4698 open (ithep,file=thetname,status='old',action='read')
4699 call getenv_loc('ROTPAR',rotname)
4700 open (irotam,file=rotname,status='old',action='read')
4701 call getenv_loc('TORPAR',torname)
4702 open (itorp,file=torname,status='old',action='read')
4703 call getenv_loc('TORDPAR',tordname)
4704 open (itordp,file=tordname,status='old',action='read')
4705 call getenv_loc('SCCORPAR',sccorname)
4706 open (isccor,file=sccorname,status='old',action='read')
4708 call getenv_loc('THETPARPDB',thetname_pdb)
4709 print *,"thetname_pdb ",thetname_pdb
4710 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
4711 print *,ithep_pdb," opened"
4713 call getenv_loc('FOURIER',fouriername)
4714 open (ifourier,file=fouriername,status='old',readonly)
4715 call getenv_loc('ELEPAR',elename)
4716 open (ielep,file=elename,status='old',readonly)
4717 call getenv_loc('SIDEPAR',sidename)
4718 open (isidep,file=sidename,status='old',readonly)
4720 call getenv_loc('THETPAR_NUCL',thetname_nucl)
4721 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
4722 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
4723 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
4724 call getenv_loc('TORPAR_NUCL',torname_nucl)
4725 open (itorp_nucl,file=torname_nucl,status='old',action='read')
4726 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
4727 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
4728 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
4729 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
4730 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
4731 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
4732 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
4733 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
4734 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
4735 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
4736 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
4737 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
4740 call getenv_loc('LIPTRANPAR',liptranname)
4741 open (iliptranpar,file=liptranname,status='old',action='read')
4742 call getenv_loc('TUBEPAR',tubename)
4743 open (itube,file=tubename,status='old',action='read')
4744 call getenv_loc('IONPAR',ionname)
4745 open (iion,file=ionname,status='old',action='read')
4748 call getenv_loc('ROTPARPDB',rotname_pdb)
4749 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
4752 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
4753 #if defined(WINIFL) || defined(WINPGI)
4754 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
4755 #elif (defined CRAY) || (defined AIX)
4756 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4758 open (iscpp_nucl,file=scpname_nucl,status='old')
4760 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
4765 ! 8/9/01 In the newest version SCp interaction constants are read from a file
4766 ! Use -DOLDSCP to use hard-coded constants instead.
4768 call getenv_loc('SCPPAR',scpname)
4769 #if defined(WINIFL) || defined(WINPGI)
4770 open (iscpp,file=scpname,status='old',readonly,shared)
4771 #elif (defined CRAY) || (defined AIX)
4772 open (iscpp,file=scpname,status='old',action='read')
4774 open (iscpp,file=scpname,status='old')
4776 open (iscpp,file=scpname,status='old',action='read')
4779 call getenv_loc('PATTERN',patname)
4780 #if defined(WINIFL) || defined(WINPGI)
4781 open (icbase,file=patname,status='old',readonly,shared)
4782 #elif (defined CRAY) || (defined AIX)
4783 open (icbase,file=patname,status='old',action='read')
4785 open (icbase,file=patname,status='old')
4787 open (icbase,file=patname,status='old',action='read')
4790 ! Open output file only for CG processes
4791 ! print *,"Processor",myrank," fg_rank",fg_rank
4792 if (fg_rank.eq.0) then
4794 if (nodes.eq.1) then
4797 npos = dlog10(dfloat(nodes-1))+1
4799 if (npos.lt.3) npos=3
4800 write (liczba,'(i1)') npos
4801 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
4803 write (liczba,form) me
4804 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
4805 liczba(:ilen(liczba))
4806 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4808 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
4810 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
4811 liczba(:ilen(liczba))//'.mol2'
4812 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4813 liczba(:ilen(liczba))//'.stat'
4815 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
4816 //liczba(:ilen(liczba))//'.stat')
4817 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
4820 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
4821 liczba(:ilen(liczba))//'.const'
4826 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
4827 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
4828 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
4829 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
4830 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
4832 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
4834 rest2name=prefix(:ilen(prefix))//'.rst'
4836 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
4839 #if defined(AIX) || defined(PGI)
4840 if (me.eq.king .or. .not. out1file) &
4841 open(iout,file=outname,status='unknown')
4843 if (fg_rank.gt.0) then
4844 write (liczba,'(i3.3)') myrank/nfgtasks
4845 write (ll,'(bz,i3.3)') fg_rank
4846 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4851 open(igeom,file=intname,status='unknown',position='append')
4852 open(ipdb,file=pdbname,status='unknown')
4853 open(imol2,file=mol2name,status='unknown')
4854 open(istat,file=statname,status='unknown',position='append')
4856 !1out open(iout,file=outname,status='unknown')
4859 if (me.eq.king .or. .not.out1file) &
4860 open(iout,file=outname,status='unknown')
4862 if (fg_rank.gt.0) then
4863 write (liczba,'(i3.3)') myrank/nfgtasks
4864 write (ll,'(bz,i3.3)') fg_rank
4865 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
4870 open(igeom,file=intname,status='unknown',access='append')
4871 open(ipdb,file=pdbname,status='unknown')
4872 open(imol2,file=mol2name,status='unknown')
4873 open(istat,file=statname,status='unknown',access='append')
4875 !1out open(iout,file=outname,status='unknown')
4878 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
4879 csa_seed=prefix(:lenpre)//'.CSA.seed'
4880 csa_history=prefix(:lenpre)//'.CSA.history'
4881 csa_bank=prefix(:lenpre)//'.CSA.bank'
4882 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
4883 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
4884 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
4885 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
4886 csa_int=prefix(:lenpre)//'.int'
4887 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
4888 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
4889 csa_in=prefix(:lenpre)//'.CSA.in'
4890 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
4893 write (iout,'(80(1h-))')
4894 write (iout,'(30x,a)') "FILE ASSIGNMENT"
4895 write (iout,'(80(1h-))')
4896 write (iout,*) "Input file : ",&
4897 pref_orig(:ilen(pref_orig))//'.inp'
4898 write (iout,*) "Output file : ",&
4899 outname(:ilen(outname))
4901 write (iout,*) "Sidechain potential file : ",&
4902 sidename(:ilen(sidename))
4904 write (iout,*) "SCp potential file : ",&
4905 scpname(:ilen(scpname))
4907 write (iout,*) "Electrostatic potential file : ",&
4908 elename(:ilen(elename))
4909 write (iout,*) "Cumulant coefficient file : ",&
4910 fouriername(:ilen(fouriername))
4911 write (iout,*) "Torsional parameter file : ",&
4912 torname(:ilen(torname))
4913 write (iout,*) "Double torsional parameter file : ",&
4914 tordname(:ilen(tordname))
4915 write (iout,*) "SCCOR parameter file : ",&
4916 sccorname(:ilen(sccorname))
4917 write (iout,*) "Bond & inertia constant file : ",&
4918 bondname(:ilen(bondname))
4919 write (iout,*) "Bending parameter file : ",&
4920 thetname(:ilen(thetname))
4921 write (iout,*) "Rotamer parameter file : ",&
4922 rotname(:ilen(rotname))
4925 write (iout,*) "Thetpdb parameter file : ",&
4926 thetname_pdb(:ilen(thetname_pdb))
4929 write (iout,*) "Threading database : ",&
4930 patname(:ilen(patname))
4932 write (iout,*)" DIRTMP : ",&
4934 write (iout,'(80(1h-))')
4937 end subroutine openunits
4938 !-----------------------------------------------------------------------------
4941 use geometry_data, only: nres,dc
4943 ! implicit real*8 (a-h,o-z)
4944 ! include 'DIMENSIONS'
4945 ! include 'COMMON.CHAIN'
4946 ! include 'COMMON.IOUNITS'
4947 ! include 'COMMON.MD'
4950 ! real(kind=8) :: var,ene
4952 open(irest2,file=rest2name,status='unknown')
4953 read(irest2,*) totT,EK,potE,totE,t_bath
4956 ! AL 4/17/17: Now reading d_t(0,:) too
4958 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
4961 ! AL 4/17/17: Now reading d_c(0,:) too
4963 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
4966 read (irest2,*) iset
4970 end subroutine readrst
4971 !-----------------------------------------------------------------------------
4972 subroutine read_fragments
4976 use control_data, only:out1file
4979 ! implicit real*8 (a-h,o-z)
4980 ! include 'DIMENSIONS'
4984 ! include 'COMMON.SETUP'
4985 ! include 'COMMON.CHAIN'
4986 ! include 'COMMON.IOUNITS'
4987 ! include 'COMMON.MD'
4988 ! include 'COMMON.CONTROL'
4991 ! real(kind=8) :: var,ene
4993 read(inp,*) nset,nfrag,npair,nfrag_back
4995 !el from module energy
4996 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
4997 if(.not.allocated(wfrag_back)) then
4998 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
4999 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5001 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5002 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5004 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5005 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5008 if(me.eq.king.or..not.out1file) &
5009 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5010 " nfrag_back",nfrag_back
5012 read(inp,*) mset(iset)
5014 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5016 if(me.eq.king.or..not.out1file) &
5017 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5018 ifrag(2,i,iset), qinfrag(i,iset)
5021 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5023 if(me.eq.king.or..not.out1file) &
5024 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5025 ipair(2,i,iset), qinpair(i,iset)
5028 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5029 wfrag_back(3,i,iset),&
5030 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5031 if(me.eq.king.or..not.out1file) &
5032 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5033 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5037 end subroutine read_fragments
5038 !-----------------------------------------------------------------------------
5040 !-----------------------------------------------------------------------------
5044 ! implicit real*8 (a-h,o-z)
5045 ! include 'DIMENSIONS'
5046 ! include 'COMMON.CSA'
5047 ! include 'COMMON.BANK'
5048 ! include 'COMMON.IOUNITS'
5050 ! integer :: ntf,ik,iw_pdb
5051 ! real(kind=8) :: var,ene
5053 open(icsa_in,file=csa_in,status="old",err=100)
5054 read(icsa_in,*) nconf
5055 read(icsa_in,*) jstart,jend
5056 read(icsa_in,*) nstmax
5057 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5058 read(icsa_in,*) nran0,nran1,irr
5059 read(icsa_in,*) nseed
5060 read(icsa_in,*) ntotal,cut1,cut2
5061 read(icsa_in,*) estop
5062 read(icsa_in,*) icmax,irestart
5063 read(icsa_in,*) ntbankm,dele,difcut
5064 read(icsa_in,*) iref,rmscut,pnccut
5065 read(icsa_in,*) ndiff
5072 end subroutine csa_read
5073 !-----------------------------------------------------------------------------
5074 subroutine initial_write
5077 ! implicit real*8 (a-h,o-z)
5078 ! include 'DIMENSIONS'
5079 ! include 'COMMON.CSA'
5080 ! include 'COMMON.BANK'
5081 ! include 'COMMON.IOUNITS'
5083 ! integer :: ntf,ik,iw_pdb
5084 ! real(kind=8) :: var,ene
5086 open(icsa_seed,file=csa_seed,status="unknown")
5087 write(icsa_seed,*) "seed"
5089 #if defined(AIX) || defined(PGI)
5090 open(icsa_history,file=csa_history,status="unknown",&
5093 open(icsa_history,file=csa_history,status="unknown",&
5096 write(icsa_history,*) nconf
5097 write(icsa_history,*) jstart,jend
5098 write(icsa_history,*) nstmax
5099 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5100 write(icsa_history,*) nran0,nran1,irr
5101 write(icsa_history,*) nseed
5102 write(icsa_history,*) ntotal,cut1,cut2
5103 write(icsa_history,*) estop
5104 write(icsa_history,*) icmax,irestart
5105 write(icsa_history,*) ntbankm,dele,difcut
5106 write(icsa_history,*) iref,rmscut,pnccut
5107 write(icsa_history,*) ndiff
5109 write(icsa_history,*)
5112 open(icsa_bank1,file=csa_bank1,status="unknown")
5113 write(icsa_bank1,*) 0
5117 end subroutine initial_write
5118 !-----------------------------------------------------------------------------
5119 subroutine restart_write
5122 ! implicit real*8 (a-h,o-z)
5123 ! include 'DIMENSIONS'
5124 ! include 'COMMON.IOUNITS'
5125 ! include 'COMMON.CSA'
5126 ! include 'COMMON.BANK'
5128 ! integer :: ntf,ik,iw_pdb
5129 ! real(kind=8) :: var,ene
5131 #if defined(AIX) || defined(PGI)
5132 open(icsa_history,file=csa_history,position="append")
5134 open(icsa_history,file=csa_history,access="append")
5136 write(icsa_history,*)
5137 write(icsa_history,*) "This is restart"
5138 write(icsa_history,*)
5139 write(icsa_history,*) nconf
5140 write(icsa_history,*) jstart,jend
5141 write(icsa_history,*) nstmax
5142 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5143 write(icsa_history,*) nran0,nran1,irr
5144 write(icsa_history,*) nseed
5145 write(icsa_history,*) ntotal,cut1,cut2
5146 write(icsa_history,*) estop
5147 write(icsa_history,*) icmax,irestart
5148 write(icsa_history,*) ntbankm,dele,difcut
5149 write(icsa_history,*) iref,rmscut,pnccut
5150 write(icsa_history,*) ndiff
5151 write(icsa_history,*)
5152 write(icsa_history,*) "irestart is: ", irestart
5154 write(icsa_history,*)
5158 end subroutine restart_write
5159 !-----------------------------------------------------------------------------
5161 !-----------------------------------------------------------------------------
5162 subroutine write_pdb(npdb,titelloc,ee)
5164 ! implicit real*8 (a-h,o-z)
5165 ! include 'DIMENSIONS'
5166 ! include 'COMMON.IOUNITS'
5167 character(len=50) :: titelloc1
5168 character*(*) :: titelloc
5169 character(len=3) :: zahl
5170 character(len=5) :: liczba5
5172 integer :: npdb !,ilen
5176 ! real(kind=8) :: var,ene
5180 if (npdb.lt.1000) then
5181 call numstr(npdb,zahl)
5182 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5184 if (npdb.lt.10000) then
5185 write(liczba5,'(i1,i4)') 0,npdb
5187 write(liczba5,'(i5)') npdb
5189 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5191 call pdbout(ee,titelloc1,ipdb)
5194 end subroutine write_pdb
5195 !-----------------------------------------------------------------------------
5197 !-----------------------------------------------------------------------------
5198 subroutine write_thread_summary
5199 ! Thread the sequence through a database of known structures
5200 use control_data, only: refstr
5202 use energy_data, only: n_ene_comp
5204 ! implicit real*8 (a-h,o-z)
5205 ! include 'DIMENSIONS'
5207 use MPI_data !include 'COMMON.INFO'
5210 ! include 'COMMON.CONTROL'
5211 ! include 'COMMON.CHAIN'
5212 ! include 'COMMON.DBASE'
5213 ! include 'COMMON.INTERACT'
5214 ! include 'COMMON.VAR'
5215 ! include 'COMMON.THREAD'
5216 ! include 'COMMON.FFIELD'
5217 ! include 'COMMON.SBRIDGE'
5218 ! include 'COMMON.HEADER'
5219 ! include 'COMMON.NAMES'
5220 ! include 'COMMON.IOUNITS'
5221 ! include 'COMMON.TIME1'
5223 integer,dimension(maxthread) :: ip
5224 real(kind=8),dimension(0:n_ene) :: energia
5226 integer :: i,j,ii,jj,ipj,ik,kk,ist
5227 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5229 write (iout,'(30x,a/)') &
5230 ' *********** Summary threading statistics ************'
5231 write (iout,'(a)') 'Initial energies:'
5232 write (iout,'(a4,2x,a12,14a14,3a8)') &
5233 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5234 'RMSnat','NatCONT','NNCONT','RMS'
5235 ! Energy sort patterns
5240 enet=ener(n_ene-1,ip(i))
5243 if (ener(n_ene-1,ip(j)).lt.enet) then
5245 enet=ener(n_ene-1,ip(j))
5257 ist=nres_base(2,ii)+ipatt(2,i)
5259 energia(i)=ener0(kk,i)
5261 etot=ener0(n_ene_comp+1,i)
5262 rmsnat=ener0(n_ene_comp+2,i)
5263 rms=ener0(n_ene_comp+3,i)
5264 frac=ener0(n_ene_comp+4,i)
5265 frac_nn=ener0(n_ene_comp+5,i)
5268 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5269 i,str_nam(ii),ist+1,&
5270 (energia(print_order(kk)),kk=1,nprint_ene),&
5271 etot,rmsnat,frac,frac_nn,rms
5273 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5274 i,str_nam(ii),ist+1,&
5275 (energia(print_order(kk)),kk=1,nprint_ene),etot
5278 write (iout,'(//a)') 'Final energies:'
5279 write (iout,'(a4,2x,a12,17a14,3a8)') &
5280 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5281 'RMSnat','NatCONT','NNCONT','RMS'
5285 ist=nres_base(2,ii)+ipatt(2,i)
5287 energia(kk)=ener(kk,ik)
5289 etot=ener(n_ene_comp+1,i)
5290 rmsnat=ener(n_ene_comp+2,i)
5291 rms=ener(n_ene_comp+3,i)
5292 frac=ener(n_ene_comp+4,i)
5293 frac_nn=ener(n_ene_comp+5,i)
5294 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5295 i,str_nam(ii),ist+1,&
5296 (energia(print_order(kk)),kk=1,nprint_ene),&
5297 etot,rmsnat,frac,frac_nn,rms
5299 write (iout,'(/a/)') 'IEXAM array:'
5300 write (iout,'(i5)') nexcl
5302 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5304 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5305 'Max. time for threading step ',max_time_for_thread,&
5306 'Average time for threading step: ',ave_time_for_thread
5308 end subroutine write_thread_summary
5309 !-----------------------------------------------------------------------------
5310 subroutine write_stat_thread(ithread,ipattern,ist)
5312 use energy_data, only: n_ene_comp
5314 ! implicit real*8 (a-h,o-z)
5315 ! include "DIMENSIONS"
5316 ! include "COMMON.CONTROL"
5317 ! include "COMMON.IOUNITS"
5318 ! include "COMMON.THREAD"
5319 ! include "COMMON.FFIELD"
5320 ! include "COMMON.DBASE"
5321 ! include "COMMON.NAMES"
5322 real(kind=8),dimension(0:n_ene) :: energia
5324 integer :: ithread,ipattern,ist,i
5325 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5327 #if defined(AIX) || defined(PGI)
5328 open(istat,file=statname,position='append')
5330 open(istat,file=statname,access='append')
5333 energia(i)=ener(i,ithread)
5335 etot=ener(n_ene_comp+1,ithread)
5336 rmsnat=ener(n_ene_comp+2,ithread)
5337 rms=ener(n_ene_comp+3,ithread)
5338 frac=ener(n_ene_comp+4,ithread)
5339 frac_nn=ener(n_ene_comp+5,ithread)
5340 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5341 ithread,str_nam(ipattern),ist+1,&
5342 (energia(print_order(i)),i=1,nprint_ene),&
5343 etot,rmsnat,frac,frac_nn,rms
5346 end subroutine write_stat_thread
5347 !-----------------------------------------------------------------------------
5349 !-----------------------------------------------------------------------------
5350 end module io_config